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ABSTRACT 

Chondrules are considered to have much information on dust particles and 
processes in the solar nebula. It is naturally expected that protoplanetary disks 
observed in present star forming regions have similar dust particles and pro- 
cesses, so study of chondrule formation may provide us great information on the 
formation of the planetary systems. 

Evaporation during chondrule melting may have resulted in depletion of 
volatile elements in chondrules. However, no evidence for a large degree of 
heavy-isotope enrichment has been reported in chondrules. In order to meet this 
observed constraint, the rapid heating rate at temperatures below the silicate 
solidus is required to suppress the isotopic fractionation. 

We have developed a new shock-wave heating model taking into account the 
radiative transfer of the dust thermal continuum emission and the line emission 
of gas molecules and calculated the thermal history of chondrules. We have 
found that optically-thin shock waves for the thermal continuum emission from 
dust particles can meet the rapid heating constraint, because the dust thermal 
emission does not keep the dust particles high temperature for a long time in the 
pre-shock region and dust particles are abruptly heated by the gas drag heating 
in the post-shock region. We have also derived the upper limit of optical depth 
of the pre-shock region using the radiative diffusion approximation, above which 
the rapid heating constraint is not satisfied. It is about 1 — 10. 

Subject headings: meteors, meteoroids — shock waves — solar system: formation 
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1. INTRODUCTION 

Chondrules are millimeter-sized, once-molten, spherical-shaped grains mainly composed 
of silicate material. They are considered to have formed from chondrule precursor dust par- 
ticles about 4.5 x 10 9 yr ago in the solar nebula; they were heated and melted through flash 
heating events in the solar nebula and cooled again to solidify in a short period of time (e.g., 
Jones et al. 2000 and references therein). So they must have great information on the early 
history of our solar system. Since it is naturally expected that protoplanetary disks around 
young stars in star forming regions have similar dust particles and processes, the study of 
chondrule formation may provide us much information on the planetary system formation. 
Chondrules have many features: physical properties (sizes, shapes, densities, degree of mag- 
netization, etc.), chemical properties (elemental abundances, degree of oxidation/reduction, 
etc.), isotopic compositions (oxygen, nitrogen, rare gases, etc.), mineralogical and petrologic 
features (structures, crystals, degrees of alteration, relicts, etc.), and so forth, each of which 
should be a clue that helps us to reveal their own formation process and that of the plane- 
tary system. Shock-wave heating model for chondrule formation has been studied by many 
authors (Hood & Horanyi 1991, 1993, Ruzmaikina & Ip 1994, Tanaka et al. 1998, Iida et 
al. 2001 (hereafter INSN), Desch & Connolly 2002 (hereafter DC02), Miura et al. 2002, 
Ciesla & Hood 2002 (hereafter CH02), Miura & Nakamoto 2005 (hereafter MN05)), and is 
considered to be one of the plausible mechanisms for chondrule formation (Boss 1996, Jones 
et al. 2000). 

Tachibana & Huss (2005) found that the degree of isotopic fractionation of sulfur in 
troilites (FeS) in chondrules from Bishunpur (LL3.1) and Semarkona (LL3.0) is quite small 
(< 0.1%/amu for all the grains, except one grain that has 0.27 ± 0.14%/amu). The absence 
of isotopic fractionation in troilites suggests that chondrules have to be heated rapidly (> 
10 4 K/hr) at a temperature range of 1273 — 1473 K, in which the isotopic fractionation should 
occur associated with evaporation of troilites (Tachibana & Huss 2005). Below 1273 K, 
troilites is solid state and it is assumed that no isotopic fractionation occurs when solid 
troilite evaporates because the time scale of sulfur diffusion in FeS is much larger than 
that of evaporation (the evaporation Peclet number for troilites at temperatures close to the 
eutectic point is about 100 for 50 /im-sized grain). On the contrary, above 1473 K, the melting 
troilite grains would have been surrounded by melted silicate. In this case, evaporation of 
sulfur would be controlled by diffusion through the surrounding silicate melt. However, since 
sulfur can hardly dissolve to the chondrule silicate melt, the troilite would not evaporate and a 
large degree of isotopic fractionation is not expected. Dust aggregates that have never melted 
before are thought to be fluffy. When such aggregates are heated, the isotopic fractionation 
of troilites should occur associated with the evaporation from inside of fluffy aggregates. 
If the heating rate is slow, it is impossible to suppress the isotopic fractionation, because 
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the duration of the evaporation becomes long enough to produce certain amount of isotopic 
fractionation. On the other hand, since once molten dust particles are not fluffy, they do 
not have to be heated rapidly to prevent from isotopic fractionation because FeS in the dust 
particle is completely covered by silicate components. Thus, we do not have to take into 
account the isotopic fractionation during the cooling phase. Namely, chondrules have to be 
heated rapidly at least in the first melting heating event. 

The isotopic fractionation can be suppressed due to not only the rapid heating but also 
the presence of back reaction from evaporated sulfur gas if the dust-to-gas mass ratio is large 
enough. However, in our situation, this effect can be negligible. Tachibana & Huss (2005) 
also calculated the degree of isotopic fractionation of sulfur under conditions of sufficiently 
high dust-to-gas mass ratio in a closed system, and concluded that the required dust-to-gas 
mass ratio to suppress the isotopic fractionation by the back reaction is higher than about 
ten thousands times the solar value. It means that the dust-to-gas mass ratio should be 
greater than about 100. It has not been well investigated whether the shock-wave heating 
(the gas drag heating) works well in such extremely high dust-to-gas mass ratio environment 
or not, and to answer this problem is beyond the scope of this paper. Thus, we do not take 
into account the back reaction in this study. 

Tachibana et al. (2004) investigated the heating rate of chondrules in the framework 
of the shock-wave heating model and concluded that the gas drag heating in the post-shock 
region can heat chondrules rapidly enough to suppress the isotopic fractionation. However, 
it is known that chondrules are also heated in the pre-shock region due to the radiation 
emitted by ambient dust particles (DC02, CH02). The effect is well known as the blanket 
effect, which was not taken into account in the study by Tachibana et al. (2004). Results 
by DC02, in which the dust thermal radiation is taken into consideration as the radiation 
source, showed that the heating speed due to the radiation is too slow to suppress the isotopic 
fractionation (~ 300K/hr). CH02 also performed numerical simulation taking the transfer 
of the dust thermal continuum emission into consideration. In contrast with DC02, results 
of CH02 showed that the heating rate of chondrules is large enough to suppress the isotopic 
fractionation even if the pre-shock region is dusty environment (dust-to-gas mass ratio is 
about 1.5, which is a few hundreds times larger than that of the solar abundance). 

In previous studies of shock-wave heating model taking into account the radiation trans- 
fer (DC02 and CH02), only the dust thermal continuum was taken into consideration as the 
radiation source. However, there is the other radiation source, the line emission of gas 
molecules. DC02 and CH02 neglected the line cooling because they assumed that the shock 
region is optically too thick to its own line emission to lose gas thermal energy effectively. On 
the contrary, MN05 showed an estimation that the post-shock gas is not so optically thick 
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and the post-shock gas in a few 100 — 1000 km behind the shock front can cool. Therefore, 
the line cooling should be taken into account. Moreover, numerical simulations in DC02 and 
CH02 have done for the limited shock conditions (the pre-shock gas number density is no — 
a few 10 14 cm" 3 and the shock velocity is v s ~ 7 — 9kms _1 ). These shock conditions are 
classified into the high-density and low-velocity shock waves in the appropriate shock con- 
dition for chondrule formation (INSN). Recently, a new shock wave generation mechanism 
which are sufficient to account for melting of dust particles was proposed (Nakamoto et al. 
2005). In the case, X-ray flares generated at close to the protosun induce the shock waves 
in the upper solar nebula where the gas density is low. The typical shock condition of the 
new mechanism is estimated as n ~ 10 11 cm -3 and v s ~ 40 — 60kms _1 . There is no study 
for the low-density and high-velocity shock waves taking into account the thermal radiation 
from dust particles. 

The purpose of this paper is to develop a new shock-wave heating code taking into 
account the radiation transfer of both the line emission of gas molecules and the dust thermal 
radiation and to investigate how the thermal history of dust particles in the shock waves is 
affected by the optical properties of the pre-shock region. Especially, we focus the heating 
rate of precursor dust particles in a temperature range of 1273 — 1473 K, in which the isotopic 
fractionation should occur. The optical properties of the flow depend mainly on the dust 
size distribution and the dust-to-gas mass ratio. Moreover, we perform the simulation with 
various shock conditions (the pre-shock gas number density and the shock velocity) in order 
to discuss which shock generating mechanism is appropriate for chondrule formation. We 
describe details of our model and basic equations in section 2. We show the calculation 
results in section 3 and discussion in section 4. Finally, we summarize our study in section 
5. 

2. MODEL AND BASIC EQUATIONS 

2.1. Overview of Model 

Basic mechanism of the shock wave heating is rather simple. Let us suppose there is a 
gas medium containing dust particles with a dynamical equilibrium, i.e. they do not have 
a relative velocity initially. And let us suppose a shock wave passes the medium. Then, 
the gas is accelerated by the gas pressure and obtains some amount of velocity, while dust 
particles tend to remain the initial position. This causes the relative velocity between the 
dust particles and the gas. When the relative velocity is present, the frictional force and 
drag heating work on the dust particles; the intensity of the force and the heating depend 
mainly on the relative velocity and the gas density. Also, the high temperature gas in the 
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post-shock region, heated by the compressional heating, heats the dust particles by thermal 
collisions. Moreover, dust particles are heated by radiation emitted from gas molecules and 
other dust particles. Dust particles are heated by those three processes, and cooled by 
emission of thermal radiation, the collision with cooler gas, and the latent heat cooling due 
to the evaporation. 

We have developed a new numerical model that simulates 1-dimensional plane-parallel 
steady shock flow including dust particles. The new model is based on our previous work 
(MN05) and added the thermal radiation transfer part to the model of MN05. We use the 
spatial coordinate, which is perpendicular to the shock front, in the shock front rest frame 
(see Fig. 1). The spatial extent of the calculation region is assumed to be the same as DC02 
for comparison (x m = 10 5 km). As will be shown later (sec. 4.3), the choice of x m is not 
a fundamental problem in this study. The origin is situated at the shock front. Initially, 
gas and dust particles flow together at the same speed (shock velocity) from x = —x m . In 
the pre-shock region, dust particles are gradually heated by the radiation as they approach 
the shock front. We focus the heating rate of the dust particles in this phase. In order 
to calculate the radiation intensity in the pre-shock region, we also have to evaluate the 
post-shock structures because the radiation mainly comes from the gas and dust particles 
in that region. When the gas flow arrives at the shock front, it is suddenly compressed 
and decelerated by the post-shock gas pressure. We consider the J-type shock and the jump 
condition is given by the Rankine-Hugoniot relation. The temperature and the density of the 
gas in the post-shock region are determined by the atomic and molecular cooling processes, 
chemical heating/cooling, and energy transfer between the gas and the dust particles. After 
passing through the shock front, dust particles have a relative velocity to the gas flow. The 
gas drag heating due to the relative velocity increases the dust temperature. In order to form 
chondrules, dust particles must be heated as high as the melting temperature. Appropriate 
shock condition for chondrule formation, in which precursor dust particles are heated enough 
to melt and do not vaporize completely, is numerically and analytically derived by INSN as 
a function of the shock velocity, v s , and the pre-shock gas number density, n . In this study, 
we simulate with twelve different shock conditions in which chondrules can be formed (see 
Fig. 2). 

Since line emission of gas molecules also heats dust particles, we have to estimate how 
much energy of line emission comes from the post-shock region. It is expected that some 
fraction of the line emission emitted from the post-shock region is absorbed by gas molecules 
before escaping toward the pre-shock region. Previous studies of shock-wave heating are 
classified into two extreme cases; one is that the line cooling is perfectly neglected (e.g., 
DC02 and CH02) and the other is that a constant fraction of line emission escapes from 
the post-shock region regardless of the column density of gas molecules (e.g., INSN). DC02 
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suggested that the shock region is optically too thick to the line emission to lose more than 
~ 10% of its energy. On the other hand, INSN considered that since the dust-melting region 
is very near from the shock front, the line cooling takes place at the region because the optical 
depth between the emitting region and the shock front is not so large. The real situation 
should be somewhere between these two models. In this study, we calculate the cooling rate 
due to the line emission taking the column density of gas molecules into consideration. Since 
a half of emitted photons goes toward upstream and the rest photons go toward downstream, 
the column density for each direction is different. Details are described in Appendix B.l. 

The gas cooling due to the line emission contributes to a part of the source term of the 
radiation field. DC02 and CH02 solved the radiation transfer of the dust thermal continuum 
emission, however, they did not take into account the line emission of gas molecules. We 
add the source term due to the line emission to the radiation model. We consider that only 
the photons which can finally escape from the post-shock region contribute the radiation 
field. This is equivalent to the on-the-spot approximation (Osterbrock 1989) used in radi- 
ation transfer and photoionization calculations for H II regions. Details of formulation is 
described in Appendix B.2. The radiation affects the thermal histories of dust particles, and 
simultaneously depends on it. Therefore, we have to calculate the gas/dust dynamics and 
the radiation transfer consistently. The procedure is iterated until the mean intensity does 
not change more than 1%. 



2.2. Gas Dynamics 

The heating and cooling processes violate the adiabatic condition of the gas. Then, 1- 
dimensional, steady flow is governed by the conservation equations of mass and momentum 
and the energy equation as follows: 

p v s = pv, (1) 

Pov* +Po = pv 2 + p, (2) 
de _ e + p dp ^ 
dt p dt 

where v s is the shock velocity, p, v, p, and e are the density, velocity, pressure of the gas and 
the gas internal energy per unit volume, respectively. The subscript "0" stands for values at 
the upstream boundary of the calculation (x = —x m ). The net cooling rate A contains the 
atomic/molecular cooling processes, chemical heating/cooling processes, and energy transfer 
between the gas and dust particles, and is given by 

A = A Lya + A H2 0(V) + AcO(V) + Ah 2 0(R) + AcO(R) + Adust + A H2 diss — r H2 form, (4) 
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where Lya, H 2 0(V), CO(V), H 2 0(R), and CO(R) stand for the cooling processes due to the 
line emission of Lyman a, H 2 vibration, CO vibration, H 2 rotation, and CO rotation, 
respectively. The energy transfer rate from gas to dust particles is A^ust- The cooling 
(heating) processes associated with the H 2 dissociation (formation) is AH 2 diss (reform)- The 
system of rate equations which governs nonequilibrium chemical reactions among gas species 
is written as 

^ 28 28 28 28 28 

= n H Y Y k ^yjVk + ^H Y Y Y k imnyiy m yn, (5) 

j=l k=l 1=1 m=l n=l 

where nn is the total number density of hydrogen nuclei, yi is the relative abundance of 
species % to nn , kjk and k\ mn are the reaction rate coefficients, respectively The first and 
second terms in Eq.(5) describe the two-body and three-body interactions in the gas phase, 
respectively. We include 156 reactions among 28 gas species. The interaction of particles 
with radiation is not included. Included processes in Eq.(5) and adopted rate coefficients 
are listed in INSN 1 . The time derivatives in Eq.(3) and (5) are Lagrangian derivatives. We 
numerically integrate Eqs.(l)-(3) and (5) using a finite difference algorithm described in 
Shapiro and Kang (1987). 

When the gas flow meets the shock front, it is suddenly compressed and decelerated by 
the pressure of gas in the post-shock region. The jump condition is given by the Rankine- 
Hugoniot relation, 

p2 vi (7 + ^)M 2 



Pl v 2 ( 7 -l)7W 2 + 2' 
T 2 _ {2 1 M 2 - (7 - 1)}{(7 - 1)M 2 + 2} 



(6) 
(7) 



T x (7 + IfM 2 

where T is the temperature of gas and 7 is the ratio of specific heat of gas. The subscript "1" 
and "2" stand for immediately before and behind the shock front, respectively. The Mach 
number in the pre-shock region is expressed as M. = V\j ^/^k^Ti/m, where fh is the mean 
molecular weight of the gas and /c B is the Boltzmann constant. 

In Eqs. (1) and (2), we neglect the mass loading by evaporation of dust materials and 
the momentum exchange with dust particles. This assumption is valid because we consider 
low dust-to-gas mass ratio (< 0.1, see sec. 2.6). However, we take into account the energy 
exchange with dust particles in Eq. (3) because the dust cooling for the gas can work well 
even in the low dust-to-gas mass ratio environment. 



1 INSN included 176 reactions among 35 gas species. In this study, we neglect D- and Si- inclusive species 
because they do not contribute to the gas cooling processes as coolant species. 
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2.3. Dust Dynamics 

If there is a relative velocity of dust particles to the gas, dust motion is changed by the 
gas drag force. The equation of motion is written by 

4 3 dv d 2 Cd I I , a s 

3 nadPd ~dT = -7ra d^-pKiKi, (8) 

where pd, Vd, C*d, and t> re i are the dust particle material density, the velocity of dust particle, 
the drag coefficient, and the relative velocity of dust particle to the gas, respectively. The 
expression of drag coefficient is given in some previous papers (e.g., Probstein 1968, Hood 
& Horanyi 1991, INSN). 

The dust temperature changes due to some heating/ cooling mechanisms; energy transfer 
with the gas, absorption of ambient radiation, radiative energy loss, and latent heat cooling 
due to the evaporation. The energy equation for the dust particle which governs the dust 
temperature is written as 

4 dT 

-iralpdC— = 47ra^(r g _ p + r rad - A rad - A cvap ), (9) 

where C and Td are the specific heat of dust particles and the dust temperature, respectively. 
Regarding the heating rate of the energy transfer with gas r g - p , we adopt the same expression 
of previous studies (Probstein 1968, Hood & Horanyi 1991). The heating/cooling rates due 
to the absorption of ambient radiation r ra( j, the emission of dust thermal continuum radiation 
A ra d, and the latent heat cooling by evaporation A evap are given as 

Trad ^abs^ ~3~ i A rac j £cmit"~SB-^d j A cvap t/cvap-^cvap; (^^) 

where e a b s (e C mit), 3 o"sb, Jevap, and L evap are the Planck- mean absorption (emission) coef- 
ficient, the mean intensity of ambient radiation, the Stefan-Boltzmann constant, the evapo- 
ration rate of dust component per unit area, and the latent heat of dust component, respec- 
tively. We assume e a b s = e em it and use the optical properties of astronomical silicate (Miyake 
& Nakagawa 1993). If the dust temperature increases as high as the boiling temperature, 
the boiling takes place in the molten droplet. After reaching the boiling point, the dust 
temperature does not change because the input energy is consumed for the phase transition 
from liquid to gas. Assuming bubbles generated in molten dust particle go away to outside 
quickly, the boiling leads to shrinkage of the dust particle (Miura et al. 2002). 

The particle radius decreases due to the evaporation from the dust surface. The evolu- 
tion of the particle radius due to evaporation is given by (Miura et al. 2002) 

^cvap / , , s 

'lit ~ ~p~d~- [ } 
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We use the evaporation rate of forsterite given by Tsuchiyama et al. (1998). 

Finally, we should mention the fragmentation of molten dust particle due to the ram 
pressure. In the post-shock region, the ram pressure acts on the surface of molten dust 
particles. Susa & Nakamoto (2002) estimated the maximum radius of molten dust particles 
above which the ram pressure dominates the surface tension and the molten dust particle 
should be fragmented. In our study, if the dust radius is larger than the maximum radius 
during its melting phase, we assume that the dust particle is divided into two molten droplets 
which have the same volume. 



2.4. Mean Intensity of Radiation Field 

We are considering the steady, one-dimensional, plane-parallel structure for the chondrule- 
forming region. The radiation transfer equation along the ray inclining with the angle 9 
against the spatial coordinate x is written as (see Eq. A4) 

^ = -(a p + a p )(l + S), (12) 

where \i is equal to cos 9, X is the frequency-initegrated specific intensity of the radiation 
field, a p (<7 P ) is the Planck-mean absorption (scattering) coefficients of dust particles, and S 
is the frequency-integrated source function. The absorption/scattering coefficients are given 
by integrating the absorption/scattering cross sections with dust radius as 

ap ) = n d 7ra 2 J Cabs )da d , (13) 

C p / J \e S cat/ 

where n d is the number density of dust particles per unit radius per unit volume and e scat is 
the planck-mean scattering efficiency The source function includes three terms; dust thermal 
continuum emission, the line emission of gas molecules, and the scattering (see Eq. A5). If 
the scattering is taken into consideration, the radiation intensity X, which is the solution of 
the radiation transfer equation, is included in the source function. Therefore, we calculate 
the radiation transfer equation iteratively until the solution converges. 

The source term of the dust thermal continuum emission j d (see Eq. A5) is given as the 
total contribution of all dust sizes; 

3d = J n d a 2 d e cmit a S BT d da d . (14) 

DC02 assumed that sub-micron sized dust particles are dynamically well coupled to the gas 
and are always in thermal equilibrium with the gas and the radiation field, however, we solve 
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the dust thermal histories of various sized dust particles (ad = 0.01 /im — 1 cm) without the 
assumption of the thermal equilibrium. 

The source term of the line emission f (see Eq. A5) is given as the summation of 
all cooling mechanisms due to the line emissions. The line emission would be absorbed or 
scattered by the gas molecules while the photon is traveling from the emitting point to a 
certain point at which we evaluate know the radiation intensity. As a result, the value of j' 
decreases as the distance between two points increases more and more. In our model, the 
absorption of the line emissions by gas molecules is considered based on Neufeld & Kaufman 
(1993) and the effect is taken into account in the term of f (see Appendix B.2). 

Here, we define the optical depth as 

/X 
(aip + <j p )dx. (15) 
■Im 

Given the incident radiation fields and the source function at all the places S, the mean 
intensity of radiation J can be found (Mihalas & Mihalas 1999): 

J(r) = ^E 2 (r) + i j( S(r') El (\r - r'\)dr' + ±fE 2 (r m - r), (16) 

where E n is the exponential integrals; E n {x) = J™ y~ n e~ xy dy, X pre and X post are given as 
the boundary conditions of radiation intensity at x — —x m and x = x m , respectively. We set 
Xpre = Xp OS t = crT 4 /n = 1.46 x 10 5 ergcm~ 2 ster -1 s _1 , which corresponds to the black body 
radiation with temperature of 300 K. Notice that we calculate the value of j g for any set of 
r and t' based on Neufeld & Kaufman (1993). We can also solve for the net flux of radiation 
energy T (Miharas & Mihalas 1999): 

T{t) = 27rX pre £ 3 (T)+27r / S(r')E 2 (r - r')dr' 

Jo 

/•Tm 

- 2tt / S{t')E 2 (t' -r)dr' - 2nl post E 3 (r m - r) . (17) 



2.5. Physical Parameters 

The initial gas temperature is assumed to be 300 K. The initial gas-phase elemental 
abundance (ratio by number to hydrogen nuclei) is taken from Finocchi et al. (1997) for the 
solar nebula: y R = 10~ 5 , y He = 9.75 x 10~ 2 , y co = 1-065 x 10~ 4 , and y H2 o = 5.978 x 10~ 4 . 
Other hydrogen nuclei are assumed to exist as hydrogen molecules. Other species = 0. 
Moreover, it is assumed that the heavy elements like Si are in the dust particles and do not 
exist in the gas phase. 
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The precursor dust particles are assumed to be composed of only forsterite with the 
mass density p mat = 3.22 gem -3 (Saxene et al. 1993). The specific heat C = 1.42 x 
10 7 ergg _1 K _1 , the latent heat of evaporation L evap = 1.12 x 10 11 ergg _1 , and the latent 
heat of boiling Lboii = 1-64 x 10 n ergg -1 (NIST WebBook are adopted). 2 For the values of 
absorption/emission/scattering coefficients (e a bs/e C mit/escat), we use the optical properties of 
astronomical silicate (Miyake & Nakagawa 1993). 

2.6. Initial Conditions for Calculations 

Appropriate shock conditions for chondrule formation, in which the precursor dust parti- 
cles are heated enough to melt and do not vaporize completely, is numerically and analytically 
derived by INSN as a function of the pre-shock gas number density no and the shock velocity 
v s . Fig. 2 shows the appropriate shock condition which is sufficient to heat the precursor 
dust particles up to 1473 — 2000 K (gray region). Based on their results, we select some set 
of the riQ and v s for the initial conditions: v s = 8, 10, 13 km s _1 for n = 10 14 cm~ 3 , t> s = 
16, 20, 25 km s" 1 for n = 10 13 cm" 3 , v s = 30, 35, 40 km s' 1 for n = 10 12 cm" 3 , and v s = 
50, 55, 60 km s _1 for uq = 10 11 cm -3 (filled circles). We abbreviate the shock condition of 
n = 10 14 cur 3 and v s = lOkms^ 1 as "nl4vl0", and so forth (see Table 1). 

The radiation field could be significantly affected by the blanket effect due to dust 
particles. Thus, the optical properties of the flow, which are mainly determined by the 
dust-to-gas mass ratio and the initial size distribution of precursor dust particles, should be 
paid attention in our study. We assume two size distribution models; power-law distribution 
and lognormal one. In the power-law one, the dust number density per unit radius per unit 
volume is given as nd(ad) oc a^" 1 , in which we assume m = 3.5 in this study. The lognormal 
one is similar to the size distribution of chondrules in ordinary chondrites (Fig. 3). We 
assume that the size range of dust particles is from O.Olum to 1cm. We divide it into 31 
bins and calculate time evolution of those particles simultaneously with gas dynamics. We 
change the dust-to-gas mass ratio in the pre-shock region as 0.01, 0.03, and 0.1. Totally, 
we have six cases as the dust particle model for each shock condition. The dust models are 
summarized in Table 2. 



2 NIST WebBook. A gateway to the data collections of the National Institute of Standards and Technology. 
Available at http://webbook.nist.gov. 
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3. RESULTS 

We first show the effect of line emission on the gas dynamics in sec. 3.1. Second, how 
the radiation field and dust thermal histories are affected by the optical depth is shown in 
sec. 3.2 and 3.3. Finally, we focus on the heating rate of dust particle in the pre-shock region 
in sec. 3.4. 



3.1. Line Cooling 

Figure 4 shows the abundance of hydrogen molecule (top), the gas temperature (middle), 
and the heating/ cooling rates for gas (bottom) in the post-shock region. The horizontal axis 
is the distance from the shock front. The shock condition is nl4vl0 in Fig. 2 and the dust 
model is L01 in Table 2. After passing through the shock front, the gas temperature jumps 
up to 4420 K. Then the dissociation of hydrogen molecule takes place and the gas cools 
rapidly. At a few tens km behind the shock front, the dissociation of hydrogen molecule 
ceases and the formation dominates the dissociation instead. Then the heating due to the 
formation of hydrogen molecule balances with the cooling due to the dissociation of it, so 
the temperature plateau is formed at T g ~ 2000 K, on which the gas temperature does not 
change significantly. In this phase, we find that the line cooling plays an important role; 
the line cooling due to the H2O vibrational emission can remove the H2 formation energy 
effectively. Thus, the abundance of hydrogen molecule increases gradually during the gas 
temperature is on the plateau. If no line cooling exists, the cooling mechanism which can 
remove the formation energy is only the dust cooling. However, its cooling rate is too low to 
cool the gas effectively. The formation energy released until all the hydrogen nuclei become 
hydrogen molecules is estimated as E = yn 2 nB_ E{ omi ~ 10 3 ergcm~ 3 , where |/h 2 ~ 0.1, 
nn ~ 10 15 cm -3 , and Ef orm = 7.17 x 10~ 12 erg is the formation energy per molecule, and 
the cooling rate by dust particles is about A^ust — 10~ 2 erg cm" 3 s _1 (see the bottom panel 
in Fig. 4), so it takes about 10 5 s to remove the formation energy. Since the gas velocity is 
about 1 kms -1 in the post-shock region, the gas would run away from the calculation region 
before their energy is removed completely. Therefore, the line emission is important for the 
gas cooling. Moreover, this result indicates that the line emission is an important radiation 
source. 
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3.2. 



Blanket Effect 



We think that it is helpful for readers to show the structure of radiation field in the 
chondrule- forming region. Because of the blanket effect, the radiation intensity strongly 
depends on the optical depth, which is specified by the size distribution and the dust-to-gas 
mass ratio C d . Fig. 5 shows the mean intensity J (top) and the radiation flux T (bottom) 
as a function of the distance from the shock front. The shock model is nl4vl0 and the dust 
models are L10 (C d = 0.10), L03 (C d = 0.03), and L01 (C d = 0.01). We find that the 
more the pre-shock region becomes optically thick, the more strongly the radiation intensity 
is amplified by the blanket effect. This result is naturally expected. The mean intensity 
immediately in front of the shock front is J = 6.36 for L10, 3.68 for L03, and 2.85 for L01 in 
unit of 10 erg cm -2 ster -1 s _1 , respectively. The temperature of dust particles which are in 
the thermal equilibrium with the radiation field is given as T ra d = {^J I osb) 1 / 4 - Hereafter, we 
call it the radiation temperature. For each case, the radiation temperature is T ra d = 1370 K 
for L10, 1195 K for L03, and 1121 K for L01, respectively. Therefore, in the case of L10, 
the dust temperature in the pre-shock region exceeds 1273 K, which is the melting point 
of FeS, and the isotopic fractionation begins to take place. Since the heating rate due 
to the blanket effect is lower than 10 3 K/hr, the isotopic fractionation will be observed if 
chondrules have been formed in such dusty shock waves. On the contrary, in cases of L03 
and L01, the radiation temperature does not exceed 1273 K in the pre-shock region, so the 
isotopic fractionation begins to take place after passing though the shock front. Behind the 
shock front, the gas drag heats dust particles very rapidly (> 10 5 K/hr) up to 1473 K, above 
which the isotopic fractionation is expected not to occur because the silicate melt surrounds 
FeS melt and suppresses its evaporation. 

The radiation flux is also an important information to understand our study. In the 
bottom panel of Fig. 5, it is found that the radiation flux is almost constant in the pre-shock 
region. In the case of L10, the pre-shock region is optically thick (r prc = 2.45), where r prc is 
defined as the optical depth of the pre-shock region; 



The reason why the radiation flux in the pre-shock region is almost constant in spite of the 
optically thick environment is the re-emission of dust thermal continuum emission (another 
explanation about the constancy of the radiation flux is shown in Appendix D). Since the 
dust temperature near the shock front is higher than that far from the shock front, the re- 
emitted dust thermal continuum emission produces net radiation flux toward upstream. This 
effect can be understood with the radiative diffusion approximation (Rybicki & Lightman 
1979). We discuss it more in sec. 4.1. 




(18) 
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3.3. Dust Temperature in Pre-Shock Region 

Figure 6 shows dust thermal histories in the pre-shock region. Calculation conditions 
are the same as those of Fig. 5; the shock condition is assumed to be nl4vl0, and dust 
models are L10, L03, and L01. A horizontal dotted line indicates 1273 K, above which the 
isotopic fractionation takes place until the dust temperature exceeds 1473 K. In cases of L03 
and L01, dust temperatures do not exceed 1273 K in the pre-shock region, so the isotopic 
fractionation does not take place in this region. After passing through the shock front, dust 
particles are rapidly heated by the gas drag heating. It takes about 1 s for dust temperature 
to increase from 1273 K to 1473 K, so the heating rates become about 10 6 K/hr in both cases. 
It is enough large to suppress the isotopic fractionation. On the contrary, in the case of L10, 
the dust temperature exceeds 1273 K in the pre-shock region at x ~ —2 x 10 4 km. This is 
due to the ambient radiation field which is amplified by the blanket effect. After that, it 
takes about 2 x 10 3 s for dust particles to meet the shock front. After passing through the 
shock front, dust temperature increases up to 1473 K within 1 s, so the duration is negligible 
compared to 2 x 10 3 s. Therefore, the heating rate between 1273 K and 1473 K is estimated 
as about 360 K/hr. It is too slow to suppress the isotopic fractionation, so it should take 
place in the case of L10. 

Here, we compare our results with that of DC02 and CH02, which took into account the 
radiation transfer of dust thermal continuum emission, but not the gas line emission (line 
cooling). The canonical case of DC02 had following parameters; the initial gas density is 
10~ 9 gcm -3 , the shock velocity is 7kms _1 , and the dust-to-gas mass ratio is 0.005. Despite 
the low dust-to-gas mass ratio, the pre-shock dust temperature reaches about 1700 K at the 
shock front (see Fig. 4 of DC02). We guess that it is because of their treatment of the 
sub-/zm sized dust particles in the source function (Eq. 5 of DC02). They used the gas 
temperature in the source term of sub-/xm sized dust particles. However, dust temperature 
is not perfectly the same as the gas temperature due to the ambient radiation field and the 
radiative cooling of dust particle themselves. For example, in the case of Fig. 4, the post- 
shock gas temperature is 2049 K at x ~ 10 3 km, where the //m-sized dust particles are in 
thermal equilibrium with the gas and the radiation field because there is no relative velocity 
between the gas and /xm-sized dust particles. The /zm-sized dust temperature at this point 
is 1247 K according to our results. The radiation source term is proportional to the forth 
power of the temperature. Thus, DC02 overestimated the source term of sub-/nn sized dust 
particles about by one order of magnitude. On the contrary, we think that the model of 
CH02 underestimated the radiation intensity in the pre-shock region because they did not 
take into account the line emission of the gas. As we showed in sec. 3.1, some amount of 
the line emission emitted in the post-shock region can penetrate into the pre-shock region. 
Considering the line emission, the radiation intensity in the pre-shock region should be larger 
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than that of the model without this effect. 



3.4. Heating Rate 

We calculated the thermal histories of dust particles for various shock conditions and 
dust models and obtained the dust temperature immediately in front of the shock front and 
the heating rate of dust particles in a temperature range of 1273 — 1473 K, in which the 
isotopic fractionation occurs. Our results are summarized in Tables 3—5. The heating rate, 
Rheat, is defined as 

1473 K- 1273 K 

-ftheat = ^ , (19) 

where At is the duration in which dust temperature increases from 1273 K to 1473 K in the 
heating phase. In Fig. 7(a), we plot the results of the dust temperature just in front of the 
shock front T sf as a function of the optical depth of the pre-shock region. In shock conditions 
nl4vl0 and nl3v20, it is found that the more the pre-shock region is optically thick, the 
higher the dust temperature just in front of the shock front becomes. If the optical depth of 
the pre-shock region r pre exceeds about 2, the pre-shock dust temperature exceeds 1273 K, 
so the isotopic fractionation takes place in the pre-shock region. On the contrary, in shock 
conditions nl2v35 and nllv55, the dust temperatures just in front of the shock front seem 
not to depend on r prc because the optical depth is very small for any dust models. The 
maximum value of r pre is 0.14 for nl2v35 and 0.013 for nllv55. In such an optically thin 
environment, the blanket effect does not work well. Fig. 7(b) shows the heating rates of 
dust particles -Rheat- It is found that if the pre-shock dust temperature exceeds 1273 K, the 
heating rate decreases drastically. In these cases, the heating rate is about a few lOOK/hr 
and it is too slow to suppress the isotopic fractionation. On the contrary, if the pre-shock 
dust temperature does not exceed 1273 K, the heating rate is very large (~ 10 6 K/hr). In 
these cases, the heating rate is determined by the gas drag heating in the post-shock region 
and the isotopic fractionation can be suppressed by the rapid heating. From these results, 
we find that the heating rates are clearly separated into two cases; the very slow heating 
(-Rheat < 10 3 K/hr) and the very rapid heating (i?h ca t ~ 10 6 K/hr). The critical optical depth 
of the pre-shock region above which the isotopic fractionation should occur is between 1.56 
and 2.45 in the shock condition of nl4vl0. 
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4. DISCUSSIONS 

In sec. 3.4, we found that the more the pre-shock region is optically thick, the higher 
the pre-shock dust temperature becomes. Moreover, we found that the heating rates of dust 
particles tend to drop drastically when the dust temperatures just in front of the shock front, 
defined as T sf , exceed 1273 K. Although this condition might not be satisfied any time, we 
think that T s { can be a good indicator in order to judge whether the rapid heating constraint 
is satisfied or not 3 . 

According to our simulation results, the critical optical depth of the pre-shock region 
above which T sf exceeds 1273 K is between 1.56 and 2.45 for shock condition of nl4vl0. In 
sec. 4.1, we derive the pre-shock dust temperature analytically using the radiative diffusion 
approximation. In sec. 4.2, we estimate the upper limit of the optical depth using discussion 
in sec. 4.1. We discuss cases of larger x m value assuming spatially extended shock waves (sec. 
4.3). We mention the cooling rate constraint for chondrule formation in sec. 4.4 and which 
shock wave generating mechanism is appropriate for chondrule formation in sec. 4.5. We 
also consider the effect of porous structures of precursor dust particles on the dust thermal 
histories in sec. 4.6. Finally, we discuss the validity of the one-dimensional plane-parallel 
geometry in sec. 4.7 and the rapid heating constraint in sec. 4.8. 



4.1. Pre-shock Dust Temperature 

In the bottom panel of Fig. 5, it is found that the radiation flux is almost constant in the 
pre-shock region. In the case of L10, the pre-shock region is optically thick (r prc = 2.45). In 
this case, it is expected that the radiation flux emitted from the post-shock region decreases 
by a factor of about e~ 2 - 45 ~ 0.086 until it reaches the pre-shock boundary (x = —x m ). The 
reason why the radiation flux is almost constant in spite of the optically thick environment 
is the re-emission of dust thermal continuum emission. The dust temperature near the shock 
front is higher than that far from the shock front, so the re-emitted dust thermal continuum 
emission produces net radiation flux toward upstream. In the optically thick environment, 
the equation of radiative diffusion (Rosseland approximation) gives the relation between the 
radiative flux and the temperature gradient as (Rybicki & Lightman 1979) 

F (20) 
3 dr v 1 



3 One of the exceptions is thought to be the shock waves with very small spatial scale, such as the 
planetesimal bow shock. We discuss that in sec. 4.5.2. 
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where dr = K^dx and is the Rosseland mean absorption coefficient. In this discussion, 
we assume /% = k, for simplicity. Assuming T to be constant and integrating Eq. (20) from 
x = —x m to x — 0, we obtain the dust temperature just in front of the shock front T sf as 

r^T 4 + -^-|% ei (21) 

40SB 

where T is the dust temperature at x — —x m and estimated by assuming the thermal 
equilibrium with the radiation field as T = (tvJ'/o'sb) 1 ^- At x — —x m , since the radiation 
intensity coming from outer region (x < — x m ) is weaker than that from the chondrule- 
forming region, we neglect its effect. Assuming the isotropic intensity in an angle range of 
7r/2 < 6 < tt, the mean intensity at x — —x m is given as J = 1/2 = \J r \/(2ir). Next, the 
radiation flux T can be estimated by a following consideration. The origin of the radiative 
energy is the gas kinetic energy entering into the shock front. The gas kinetic energy flux 
is approximately estimated as \p§v\ . Some fraction of the gas energy flux returns toward 
upstream as a form of radiation flux. Assuming the fraction is /, we obtain \T\ = fpo^s- 
Our calculation results of / are summarized in Tables 3-5. Finally, we can estimate the dust 
temperature immediately in front of the shock front as 

In Fig. 8, we compare the pre-shock dust temperature estimated by Eq. (22) with the 
calculated results for various shock conditions: nl4vl0 (top left), nl3v20 (top right), nl2v35 
(bottom left), and vllv55 (bottom right), respectively. We substitute po = fan® into Eq. 
(22), where the mean molecular weight m = 2.35 x lCT 24 g resulting from the initial gas 
abundance we adopted in this study and draw three curves which correspond to / = 0.3, 
0.5, and 1.0, respectively. Triangles indicate the calculation results for the log-normal dust 
size distribution and circles are for power-law one. The calculation results seem to well match 
with the analytic formula (Eq. 22) in cases of r prc > 1 if the difference of / among each 
dust model is taken into consideration. On the contrary, in the optically thin environment 
(Yprc < 1), especially in the shock conditions of nl2v35 and nllv55, T s f of calculation results 
are larger than a curve indicating / = 1 which means that all of the gas energy flux flowing 
into the shock front is converted into the radiation energy flux and it returns upstream. The 
reason why the calculation results exceed the upper limit is thought to be due to the optically 
thin environment in which the radiative diffusion approximation cannot be applied. However, 
our results show that the analytic formula Eq. (22) is valid for the relatively optically thick 
environment (r prc > 1) in which the radiative diffusion approximation can be used. 

We would like to mention the effect of the radiation intensity coming from outer region 
(x < —x m ), which is neglected in above estimation. In an optically thick environment, the 
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radiation intensity is proportional to T£. We assume the radiation intensity coming from 
outer region corresponds to the black body radiation of 300 K. If the typical dust temperature 
in chondrule-forming region is 1200 K, the ratio of the radiation intensity coming from outer 
region to that in the chondrule-forming region becomes only about 0.004. Therefore, our 
assumption that the radiation intensity coming from outer region is negligible is valid. Even 
if the temperature in outer region is 500 K, its contribution is only about 0.03 in radiation 
intensity. To summarize, the dust temperature immediately in front of the shock front T sf 
does not depend on the temperature in outer environment significantly. 



4.2. Critical Optical Depth 

In order to suppress the isotopic fractionation, it seems that the pre-shock dust temper- 
ature, T s f in Eq. (22), should be lower than 1273 K. The requirement gives the upper limit 
of the optical depth of the pre-shock region above which the isotopic fractionation should 
occur. Rewriting Eq. (22), the condition in which the isotopic fractionation is suppressed is 
given as 

4a SB (1273K) 4 2 

Tpre < Tex = g • (23) 

We plot the upper limit of the optical depth r cr for each shock condition as a function of / 
in Fig. 9. In all shock conditions, large amount of radiation flux returning from the shock 
front requires the smaller optical depth to suppress the isotopic fractionation. Moreover, 
the shock condition nllv55 allows the larger optical depth than that of nl4vl0 because the 
radiation flux is smaller than that of nl4vl0. If the optical depth of the pre-shock region is 
lower than the upper limit in Fig. 9, the pre-shock dust temperature does not exceed 1273 K 
and the isotopic fractionation would be suppressed in the pre-shock region. In this case, in 
the post-shock region, since the silicate component of dust particles melts quickly due to the 
gas drag heating, the evaporation of FeS is suppressed and the isotopic fractionation does 
not take place. To summarize, the upper limit of the pre-shock optical depth r prc is about 
1 — 10 for various shock conditions (no = 10 11 — 10 14 cm~ 3 ). 



4.3. Different Spatial Scale 

In this study, we assumed x m = 10 5 km, which indicates that the spatial scale of shock- 
wave heating region is 10 5 km. In other words, only dust particles located within 10 5 km 
from the shock front contribute to the blanket effect, and the thermal radiation of dust 
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particles in outer region runs away effectively. The length of 10 5 km corresponds to about 
0.001 AU. On the contrary, if shock waves were generated by the marginally gravitational 
instability (e.g., Boss & Durisen 2005), the extent of chondrule-forming region would be more 
larger. However, as we pointed in sec. 4.1, the dust temperature just in front of the shock 
front is determined by the gas energy flux flowing into the shock front and the pre-shock 
optical depth r prc . Since T s f does not explicitly depend on the dimension of the pre-shock 
region, the spatial scale of the shock waves is not an important factor to consider the rapid 
heating constraint. In this section, in order to confirm that the optical depth of the pre-shock 
region r prc is a fundamental parameter to determine the pre-shock dust thermal profile, we 
numerically calculate the case of a different spatial scale (x m = 10 6 km) and the dust-to-gas 
mass ratio d = 0.01, which results into almost the same r pre as the dust-to-gas mass ratio 
Cd = 0.10 in the standard spatial scale (x m = 10 5 km). 

Figure 10 shows the dust thermal histories in the pre-shock region for different spatial 
scale: x m = 10 6 km and Cd = 0.01 (top), and x m = 10 5 km and Cd = 0.10 (bottom). 
The shock condition is nl4vl0 and the dust size distribution is power-law, which result 
into the pre-shock optical depth of r prc = 10.5 for top panel and r prc = 14.1 for bottom 
panel, respectively. Two dashed lines show the temperature range in which the isotopic 
fractionation takes plane. It is found that the dust thermal profiles in two cases are very 
similar in spite that the spatial scales are ten times different. These results indicate that 
the pre-shock optical depth r prc is an fundamental indicator to determine the dust thermal 
histories in the pre-shock region. 

Strictly speaking, these two dust thermal profiles are not the same exactly because of 
the slight difference in the pre-shock optical depth r prc . It is due to the difference of the 
duration for the dust particles to pass through the pre-shock region; in the top panel, it 
takes ten times longer than that of the bottom panel. As a result, even larger dust particles 
can evaporate completely in the top panel. In the top panel, the dust particles whose radii 
are smaller than 2.51 /xm evaporate completely in the pre-shock region, on the contrary, in 
the bottom panel, the dust particles of 0.63 fim in radius can survive in this region. As a 
result, the pre-shock optical depth in the top panel becomes smaller than that of the bottom 
panel. The difference in the optical properties in the pre-shock region causes the difference 
in the dust thermal profiles in that region. 

4.4. Cooling Rate Constraint 

For the chondrule formation, not only the rapid heating constraint, but also the appro- 
priate cooling rate should be satisfied. In our simulations, we can obtain the thermal histories 
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of the dust particles for various radii (peak temperature T pea k, heating rate i?heat, and cool- 
ing rate -R coo i) as long as the dust particles do not evaporate away in the chondrule-forming 
region. 

4-4- 1- Appropriate Cooling Rate 

In order to produce chondrules, it is believed that cooling rates after the flash melting 
event should be appropriate values. However, we discussed the appropriate shock waves 
only in the view of the rapid heating constraint in the previous section. It is because the 
cooling rate constraint has not been well determined and could not restrict the flash heating 
mechanisms so strongly. The data of Jones & Lofgren (1993), who compared the textures 
and zoning profiles of olivine grains tens of micrometers across in natural chondrules and 
experimental analogs, clearly showed that the range of cooling rates experienced by type IIA 
chondrules was 5 — 100 K/hr. For higher cooling rates (500 — 1000 K/hr), olivine textures in 
these experiments were generally more elongate and skeletal than those grown at the slower 
cooling rates, and moreover, chemical zoning was rather limited. On the contrary, Wasson 
& Rubin (2003) measured overgrowths on low-FeO relict grains contained within almost 
all of type II chondrules and found that the overgrowths are narrow, in the range of 2 — 
12 /iin. Wasson (2004) inferred the cooling rate of chondrules assuming that the overgrowths 
represent the maximum amount of material that can reach the crystal/melt interface by 
diffusion in the melt during the cooling interval. Diffusion length can be estimated as 5 2 ~ 
Dt coo \, where 5 is the diffusion length, D is the diffusion coefficient, and t coo i is the cooling 
time scale. Thus, the results that the amount of the overgrowths was about 30 times smaller 
than the grain sizes simulated in order to estimate cooling rates {e.g., Jones & Lofgren 
1993) indicate that the cooling rates should be about 900 (= 30 2 ) times faster than that 
suggested by Jones & Lofgren (1993). Additionally, Yurimoto & Wasson (2002) modeled 
the O-isotopic and FeO/(FeO+MgO) gradients in a type II chondrule from the Yamato- 
81020 CO3.0 chondrite and inferred an extremely fast cooling rate (10 5 — 10 6 K/hr) at high 
temperatures (~ 1900 K). To summarize, though the information of the cooling rates is 
very important for chondrule formation mechanism, the appropriate values have not been 
well determined yet. We would like to discuss the appropriate shock waves for chondrule 
formation in the view of the rapid heating and the appropriate cooling rate constraints in 
the future paper. 



- 21 - 



4-4-2. Heating Rate vs. Cooling Rate 

In order to discuss the chondrule formation, the cooling rate in the phase that the 
precursor dust particle re-solidifies after the melting phase is an important information. We 
numerically calculate the thermal histories of the dust particles in the shock-wave heating. 
We estimate the cooling rate as 

_ 1550 K- 1400 K 

-Kcool = , 

where At is the duration in which dust temperature decreases from 1550 K to 1400 K in the 
cooling phase. The results for various shock conditions and the dust models are summarized 
in Tables 3—5. When the dust particle evaporates completely (final dust radius = jum), 
we could not obtain the peak dust temperature T pea k, the heating rate Rheax, and the cooling 
rate i? coo i- We also could not obtain the heating rate -Rhcat if T pea ^ < 1473 K and the cooling 
rate -R coo i if ^peak < 1550 K. In those cases, data are not listed in the tables. 

Figures 11-13 show the plot of the heating rate vs. the cooling rate. Symbols indicate 
the simulation results and the different symbols represent different gas number density in 
the pre-shock region; no = 10 14 cm -3 (circle), 10 13 cm" 3 (square), 10 12 cm~ 3 (triangle), and 
10 11 cm -3 (inverse triangle), respectively. The vertical dashed line indicates the lower limit 
of the heating rate below which the isotopic fractionation would take place. We draw the 
horizontal lines at -R CO oi = 5 K hr _1 and 100 K hr _1 between which the textures and zoning 
profiles of olivine grains tens of micrometers across in natural chondrules are well explained 
(Jones & Lofgren 1993). We also draw the horizontal line at a higher cooling rate (-R CO oi > 
5000 Khr" 1 ) suggested from the measurement of the overgrowth on low-FeO relict grains 
contained within almost type II chondrules (Wasson 2004). 

It is found that the numerical results are roughly classified into two categories; the first is 
the slow heating (i?heat = 10 2 — 10 3 K hr -1 ) and slow cooling (-R CO oi = 5 — 20 K hr" 1 ), and the 
second is the rapid heating (-Rhcat ~ 10 6 K hr -1 ) and rapid cooling (-R coo i = 10 3 — 10 5 K hr" 1 ). 
The former case appears when the gas number density in the pre-shock region is relatively 
high (no > 10 13 cm -3 ) and the optical depth is larger than about unity. On the contrary, 
if the pre-shock gas number density is low (no < 10 12 cm" 3 ) or the optical depth is lower 
than about unity, the thermal histories of the dust particles become the latter case. Namely, 
the heating rate and the cooling rate become extremely small only when the radiation field 
is strong enough. The reason can be easily understood as follows. If the mean intensity 
of the radiation field is strong enough to heat the dust particle up to 1273 K or more in 
the pre-shock region, the heating rate is determined at the pre-shock region and it becomes 
small (see sec. 3.4). Moreover, when the radiation field is strong in the pre-shock region, 
the mean intensity of the radiation field is also strong in the post-shock region (see Fig. 5). 
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The strong radiation field keeps the dust particle at high temperature for a long time, so the 
cooling rate becomes small. 

The rapid heating/cooling case satisfies the rapid heating constraint for suppressing 
the isotopic fractionation. The cooling rate meets the measurement of the overgrowths on 
low-FeO relict grains contained within almost all of type II chondrules (Wasson 2004). On 
the contrary, in the slow heating/cooling case, the isotopic fractionation must take place. 
However, the cooling rate well matches with the estimation suggested by comparisons of the 
textures and zoning profiles of olivine grains tens of micrometers across in natural chondrules 
and experimental analogs (Jones & Lofgren 1993). 

4.5. Appropriate Shock Mechanism 

In order to form chondrules, the shock wave generation mechanism should satisfy at 
least the rapid heating constraint. We discuss which shock wave generation mechanism is 
appropriate for chondrule formation in view of the rapid heating constraint. 

4-5.1. Gravitational Instability 

Recently, Boss & Durisen (2005) showed that the marginally gravitational instability in 
the nebula can drive inward spiral shock fronts at asteroidal orbits, sufficient to account for 
melting of dust particles. Shock waves generated by this mechanism are similar to the shock 
condition of nl4vl0. They did not examine if the rapid heating is satisfied or not. If the 
spatial scale of shock waves is 10 5 km, which is the same as the calculation region we assumed 
in this study, a shock wave with the dust-to-gas mass ratio less than about 0.03 satisfies the 
rapid heating condition. However, if the spatial scale is 10 6 km, which is ten times larger 
than that of our assumption, it is impossible to suppress the isotopic fractionation even if 
the dust-to-gas mass ratio is 0.01, which is a standard value of the minimum mass solar 
nebula model. We clearly showed this conclusion in Fig. 10 for the dust model POL For a 
case that all dust particles have the same radius, we can easily derive the same conclusion as 
follows. The optical depth of the pre-shock region is calculated as r pre ~ (a p + o~ p )x m , where 
Oi p + Op = vra^(e a bs + e scat )nd. Here, we assume that the dust particles have the same radius 
of 500 /j,m. Assuming p = 4 x 10~ 10 gcm~ 3 and the dust-to-gas mass ratio Cd = 0.01, the 
number density of dust particles is n d ~ 2 x 10~ 9 cm~ 3 . For the astronomical silicate with 
radius of 500 /xm, e a b s + e scat ~ 2. Therefore, substituting x m = 10 6 km, we finally obtain 
Tpre ~ 3 and it exceeds the upper limit of the optical depth in Fig. 9. 
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4-5.2. Planetesimal Bow Shock 

Near the nebula midplane, possible energy sources for forming chondrules include bow 
waves upstream of planetesimals scattered gravitationally into eccentric/inclined orbits by 
proto- Jupiter (Hood 1998, Weidenschilling et al. 1998). The gravitational instability shock 
wave mechanism stated above predicts relatively large-scale sizes (> 10 4 km) for shock- 
heated regions, while the planetesimal bow shock predicts much smaller scale size (~ 10 km). 
Because of the small scale size, the optical depth of the chondrule-forming region will be very 
small (r prc ~ 1CT 5 ). Since the blanket effect does not work well in such an optically thin 
environment, the isotopic fractionation would be prevented. 

Even if the optical depth exceeds r cr , the rapid heating constraint might be satisfied in 
such small spatial scale. In order to prevent the isotopic fractionation, the precursor dust 
particles should be heated from below 1273 K to above 1473 K within ~ (1473—1273) / -Rheat = 
72 s, where we substitute -Rheat — 10 4 K hr _1 . If we take the shock velocity of v s = 10 km s _1 , 
the precursor dust particles can travel 720 km within 72 s. On the contrary, the spatial 
scale of the planetesimal bow shock is roughly estimated as 10 km, it does not take 72 s to 
pass through the pre-shock region. Therefore, the isotopic fractionation would not occur 
significantly. However, since the optically thick environment with such small spatial scale 
requires a very large dust-to-gas mass ratio (Cd = 100 or more), this situation might not be 
present in the protoplanetary disk. 

4-5.3. Clumpy Cloud Accretion 

Boss & Graham (1993) suggested that the source of the shocks possibly responsible 
for chondrule formation was episodic accretion onto the solar nebula of low-mass clumps. 
The existence of the clumps is suggested by some observational evidences (Bellingham & 
Rossano 1980). They considered that the clumps are orbiting above the protoplanetary disk 
and sometimes obscuring the central star. When the clumps accrete onto the disk surface, 
the shock waves are generated at the upper region of the disk. Boss & Graham assumed 
that the following parameters characterize a "standard" clump: mass M c = 10 22 g; density 
p c = 3 x 10~ 15 gcm~ 3 ; and radius R c = 9 x 10 11 cm. Assuming the same parameters as the 
case of the gravitational instability except for x m = 10 km and po = 4 x 10~ 13 gcm -3 (upper 
region of the protoplanetary disk is low density), we obtained the pre-shock optical depth of 
T pre ~ 0.3. The optical depth is low enough to prevent the isotopic fractionation. Therefore, 
the accretion shock waves seem to match with the rapid heating constraint. However, Tanaka 
et al. (1998) carried out the hydrodynamic simulation of the accreting clumps based on the 
standard clump and they concluded that more than about 10 2 to 10 3 times surface densities 
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for the same size of the representative clump are needed to produce chondrules. They also 
suggested that there are not so many observational studies of the clumpy clouds around 
the young stellar objects, so more massive clumpy clouds which are sufficient to melt the 
precursor dust particles might exist. 

4-5.4- X-ray Flare Induced Shock Wave 

Nakamoto et al. (2005) proposed a new shock wave generation mechanism: shock waves 
in the upper region of the solar nebula induced by X-ray flares associated with the young 
Sun. X-ray flares, common among T Tauri stars (Feigelson et al. 2002), emit plasma gas, 
which cools to be a strong neutral gas wind. Hayashi et al. (1996) suggested that the 
magnetocentrifugal acceleration (Blandford & Payne 1982) also plays an important role to 
produce the wind. Because of the enormous amount of energy released by the X-ray flares 
and the magnetocentrifugal acceleration, the flares should have some effects on the dynamics 
and energetics of a protoplanetary disk around the star. Nakamoto et al. (2005) carried 
out 2-D MHD numerical simulations of X-ray flares around a central star and expanding 
magnetic bubles/winds with a disk. They showed that shock waves which are sufficient to 
melt precursor dust particles can be generated in the upper region of the nebula inside about 
the asteroid belt. In that case, it is thought that the shock waves are generated in a very wide 
region at the upper surface of the disk. The neutral gas wind travels in almost parallel with 
the disk surface much faster than the shock waves propagating vertically in the disk toward 
the midplane. Therefore, the shock front caused by the neutral gas wind is approximately 
parallel to the disk surface and the horizontal scale of the shock front would be larger than 
the scale height of the gas disk. It indicates that the dust particles located at the midplane 
also can contribute to the blanket effect. As a result, the optical depth perpendicular to 
the shock front might become much larger than unity. To summarize, there is a possibility 
that the X-ray flare induced shock waves result in the isotopic fractionation in spite that 
the chondrule formation occurs in the low gas density region. However, if the shock wave is 
non-steady, this difficulty may be avoided, because the dust in the midplane can become a 
large absorber of the radiation flux rather than a reflector and the blanket. This possibility 
will be discussed in the following subsection. 

4-5.5. Other Possibilities 

We verified four shock generation models as the chondrule formation mechanisms. We 
found that most of shock waves can be classified into two categories; the rapid heating and 
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the rapid cooling, or the slow heating and the slow cooling (see sec. 4.4.2). If supporting an 
opinion that the slow cooling rate is appropriate for producing chondrules (Jones & Lofgren 
1993), there are few cases in which both of two constraints are satisfied simultaneously (see 
Fig. 11 — 13). How can we explain the chondrule formation with the rapid heating and the 
appropriate slow cooling? 

One possibility to avoid above difficulties is to consider the non-steady shock wave 
model. When shock waves are generated in the solar nebula, the shocked gas just behind the 
shock front emits the line photons and they penetrate into the pre-shock region. Although 
the radiation heats the dust particles in the pre-shock region, it will take time before the dust 
temperatures have a steady profile. If the shock waves disappear within a shorter period 
of time, the dust temperatures would not be as high as the steady profiles. Therefore, the 
rapid heating constraint might be satisfied even in a relatively optically thick environment. 
The cooling rate of the precursor dust particles has not been investigated in the framework 
of the non-steady shock waves, so it is an important issue which should be studied in the 
future. 

Another possibility is a multiple shock heating scenario. In this study, we assume the 
porous structures for the precursor dust particles in order to consider the evaporation of 
sulfur inside the aggregates. In that case, the rapid heating constraint is required in order 
to prevent the isotopic fractionation. However, if the precursor dust particles have been 
once melted by some heating process, since the evaporation of sulfur inside the aggregates 
can hardly occur, the rapid heating might not be required. Based on these discussions, we 
can consider the multiple shock heating scenario in which the rapid heating constraint and 
the appropriate slow cooling are both satisfied. In the first heating event, the precursor 
dust particles should be heated by the shock waves in which the rapid heating constraint 
is satisfied. At that time, the cooling rate can take any values. After that, in the final 
heating event, the once melted precursor particles should be heated by the shock waves with 
the appropriate slow cooling. At that time, the rapid heating is not required because the 
sulfur cannot evaporate away from inside of the dense particles. In fact, there are some 
observational evidences of the multiple heating in chondrule formation processes (Jones et 
al. 2000), so the above scenario seems to be reasonable. However, it is not clear whether the 
isotopic fractionation can take place from the once- molten particles or not. The condition 
of the isotopic fractionation from once-molten particles is also an interesting study. 
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4.6. Porosity of Precursor Dust Particles 

For calculating the dust dynamics, we assumed dense structures for the chondrule pre- 
cursor dust particles (mass density is the same as a value of forsterite). However, it is thought 
that the precursor dust particles before melting have porous structures, therefore sulfur can 
evaporate from inside the dust particles and the isotopic fractionation can take place. If the 
porous structures are considered, a cross section of the dust particle would be larger than 
the dust model that we adopted in this study. As a result, the optical depth of the pre-shock 
region r prc would be larger, even if the dust model (dust size distribution, dust-to-gas mass 
ratio) is fixed. In this section, we consider to what extent the porous structures have an 
effect on our conclusions. 

Figure 14 shows schematic pictures of precursor dust particles having various porous 
structures. The dark gray region indicates a dense region in which the porosity is 0%. The 
light gray region indicates a porous region in which the porosity is assumed to be 80% so that 
sulfur can evaporate from inside of the porous region. We assume that the porous region is 
surrounding around the dense core. The volume percent of each region is displayed below 
each picture. We can obtain the dust radius of each structure from the volume percent of 
each region and the porosity of the porous region, assuming that the total masses of those 
particles are the same. We also display the radius and the cross section in units of the values 
for the case (a), in which all the region of the precursor dust particle is dense. 

If all the region of the precursor dust particles is porous (case (d)), the radius and the 
cross section are about 1.7 and 2.9 times larger than that of the case (a), respectively. It 
indicates that the optical depth of the pre-shock region r prc becomes about 3 times larger 
than that obtained in this study. For example, in the case that the shock model is nl4vl0 
and dust model L03 (see Table 4), we obtained r prc = 0.73. However, if all dust particles 
have porous structures like the case (d) in Fig. 14, r pre would be about 2.13. The increase 
of T pre would result in that the dust temperature just in front of the shock front T s f exceeds 
1273 K (see Fig. 8) and the isotopic fractionation takes place. On the contrary, in a case 
of the dust model L01, the optical depth of the pre-shock region r prc would be about 0.7 
if all precursor dust particles are assumed to have porous structures like the case (d) in 
Fig. 14. This expected value of r pre seems not to be so large that the isotopic fractionation 
would not occur. To summarize, we think that the discussion about the upper limit of r prc 
for preventing the isotopic fractionation in sees. 4.1 and 4.2 is not affected significantly by 
considering the porous structures because the dust temperature just in front of the shock 
front T s f mainly depends only on r pre and the net flux of the radiation T (see Section 4.1). 
However, the corresponding dust model (dust size distribution, dust-to-gas mass ratio) for a 
certain value of r nro would be affected. 
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Finally, we notice that the structures of chondrule precursor dust particles have not been 
well investigated. Although we assumed the porosity of the porous region as 80% for sulfur 
to evaporate from inside, a value of the porosity is unknown. Moreover, it is not certain 
whether all region of the precursor dust particles should be porous or not. We can imagine 
the mixed structure like the case (b) in Fig. 14, which has the dense core surrounded by the 
porous mantle. Of course, since sulfur would be able to evaporate from inside of the porous 
mantle, it must be heated rapidly enough to prevent the isotopic fractionation. In the case 
(b), the cross section is about 1.5 times larger than that of the case (a). Therefore, although 
the internal structure of precursor dust particle would affect the optical properties of the 
pre-shock region, it is a very complex problem. If we take into account the internal structure 
of the dust particle into our model, it will require a large amount of parameter sets which 
characterize the dust structure. It is beyond the scope of this paper to investigate for all of 
the parameter sets. 



4.7. Plane-parallel Geometry 

We assumed the one-dimensional, plane-parallel, and steady shock structures in this 
study. In these situations, the radiation flux JF should be constant in the pre-shock region 
because there are not any radiative sources or sinks (see Appendix D). It is clearly shown 
in our numerical results (see Fig. 5). It indicates that the radiative heating takes place 
everywhere in the pre-shock region, even at the upstream boundary of the computational 
domain (x = —x m ). The dust temperature at the boundary can be estimated from Eq. (22) 
by substituting r prc = 0. We obtain T d | x= _ Xm = 969 K for / = 0.5, p = 10~ 10 gcm~ 3 , and 
v s = lOkms -1 . On the contrary, we assume that the ambient temperature for the gas and 
the dust particles is 300 K. How can we recognize this temperature jump? 

The answer is the three-dimensional effect of the chondrule-forming region. It is con- 
sidered that the one-dimensional plane-parallel approximation would be adoptable near the 
shock front. In this region, the radiation would travel perpendicularly to the shock front 
and the radiative flux T should be constant. However, far from the shock front, which is 
roughly estimated by the spatial extent of the shock front, the plane-parallel assumption 
breaks down. The radiation can diffuse away not only perpendicularly to the shock front 
but also to the other directions, therefore, the effect of the radiative heating becomes weak 
rapidly as it keeps away from the shock front. Finally, the dust temperature becomes as the 
same as the ambient temperature. 

If the thermal histories of the dust particles below about 1000 K are in interest, it would 
be needed to develop the two- or three-dimensional shock-wave heating model. In the case of 
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the planetesimal bow shock, some studies about the multi-dimensional shock-wave heating 
have been carried out (Hood 1998, Weidenschilling et al. 1998), in which the planetesimals 
are relatively small (~ 10 — 100 km) and the multi-dimensional effect should be considered. 
In our study, we consider relatively large scale shock waves (~ 10 5 km or more). Moreover, 
the thermal histories of the dust particles below 1000 K do not play an important role for 
our purpose. Therefore, we do not treat the thermal evolutions of the dust particles before 
they enter the one-dimensional plane-parallel region. 

4.8. How Strong is Heating Constraint? 

In this study, we discussed the appropriate shock conditions and the dust models for 
chondrule formation with the rapid heating constraint, derived from the measurements of 
the isotopic composition of sulfur (Tachibana & Huss 2005). However, as commented in 
their abstract, this constraint is based on the assumption that the troilite grains within 15 
ferromagnesian chondrules they measured are primary. They considered that some troilite 
grains have survived chondrule-forming high-temperature events because they are located 
inside metal spherules within chondrules. According to them, such an occurrence is unlikely 
to be formed by secondary sulfidization processes in the solar nebula or on parent bodies. 

If those troilite grains are not primary, the rapid heating constraint would not be re- 
quired necessarily. In our study, we found that the rapid heating (i?heat ~ 10 6 Khr _1 ) does 
not seem to be compatible with the slow cooling (-R coo i ~ 5 — 20 K hr _1 ), which is considered 
by many researchers to be appropriate for chondrule formation (e.g., Jones & Lofgren 1993). 
However, we do not mean that we should adopt the rapid heating constraint rather than 
the slow cooling rate. We just investigated the thermal evolutions of the dust particles in 
the shock-wave heating for various situations by using a more realistic physical model. The 
rapid heating constraint is not needed to be assumed in order to derive our numerical results. 
We think that more experiments and observations are required for elucidating the riddles of 
chondrule formation. We believe that our simulations give a great insight of the shock-wave 
heating model for chondrule formation to readers. 

5. CONCLUSIONS 

We numerically simulated the shock-wave heating for chondrule formation taking into 
account the radiation transfer of the line emission of gas molecules and the dust thermal 
continuum emission. Regarding the line cooling due to the gas molecules, we estimated 
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the cooling rate using the photon escape probability method in order to reflect the column 
density of gas molecules. We focused the dust thermal histories in the pre-shock region, 
especially the heating rate of chondrules -Rheat, which has to be large enough to suppress the 
isotopic fractionation. In order to clear the condition in which the rapid heating constraint 
is satisfied, we performed simulations for various shock conditions, the pre-shock gas number 
density (no = 10 11 — 10 14 cm~ 3 ) and the shock velocity (v s = 10 — 55kms _1 ), and the dust 
models, the initial dust size distribution (the power-law size distribution and the lognormal 
one) and the dust-to-gas mass ratio (C d = 0.01, 0.03, and 0.10). We found the following 
results: 

1. The line emission can remove the post-shock gas thermal energy away even if the 
pre-shock gas number density is relatively high (no = 10 14 cm -3 ). Therefore, the line 
emission plays an important role as the gas cooling mechanism. Moreover, it indicates 
that it is also an important radiation source term. 

2. When the optical depth of the pre-shock region increases, it results into the higher 
dust temperatures in the pre-shock region. If the dust temperatures just in front of 
the shock front exceed 1273 K, the heating rate of the dust particles in a temperature 
range of 1273—1473 K is too slow to prevent the isotopic fractionation. On the contrary, 
the dust temperatures just in front of the shock front are lower than 1273 K, the dust 
particles are heated by the gas frictional heating in the post-shock region so rapidly 
that the isotopic fractionation is prevented. Therefore, the condition to prevent the 
isotopic fractionation is that the dust temperatures just in front of the shock front do 
not exceed 1273 K. 

3. We analytically derived the dust temperature just in front of the shock front using the 
theory of the radiative diffusion. The analytic solution well explains the results of our 
numerical simulations for optically thick case (optical depth of the pre-shock region 

Ve £ !)■ 

4. There is the upper limit of the optical depth of the pre-shock region above which the 
isotopic fractionation will occur. The value of the upper limit depends on the shock 
condition. For the low- velocity and high-density shock waves, the upper limit is about 
unity. For the high-velocity and low-density shock waves, it is about 10. 

5. The fundamental factors to determine the pre-shock dust thermal histories are the gas 
energy flux flowing into the shock front and the optical depth of the pre-shock region. 
Even if the spacial dimension of the shock- wave heating region is different, the dust 
temperatures just in front of the shock front become almost the same as long as the 
pre-shock regions take the similar values of the optical depths. 
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6. We also evaluated the cooling rate of the precursor dust particles -R coo i at a crystal- 
lization temperature. It was found that the dust thermal histories obtained by our 
simulations are roughly classified into two categories: the first is the slow heating 
(-Rhcat — 10 2 — 10 3 Khr _1 ) and slow cooling (-R coo i = 5 — 20Khr~ 1 ), and the second is 
the rapid heating (Rheat ~ 10 6 K hr _1 ) and the rapid cooling (R coo i — 10 3 — 10 5 K hr _1 ). 
If supporting an opinion that the slow cooling rate is appropriate for producing chon- 
drules, it seems difficult to meet two conditions simultaneously in the framework of a 
one-dimensional, plane-parallel, and steady shock wave model. In order to satisfy the 
rapid heating constraint and the appropriate slow cooling constraint simultaneously, a 
non-steady shock wave model or a multiple heating scenario might be needed. 
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A. Radiative Transfer Equation 

One of the purposes of this study is to develop the numerical code of shock-wave heating 
model taking into account the dust thermal continuum emission and the line emission of gas 
molecules. However, to solve the transfer of the line emission with no approximation needs 
a huge amount of computational time. So, we approximate the line transfer problem as 
follows. 

Assuming the isotropic scattering, the frequency-dependent radiation transfer equation 
taking into account absorption and scattering by the gas molecules and dust particles is 
given as (e.g., Rybicki & Lightman 1979) 

= -{ol v , & + a vA + a Vjg + o vA )T v + j uA + j^g + (<7„ )g + a vA )J v , (Al) 

where a u , <r u , and j u are the absorption coefficient, the scattering coefficient, and the emis- 
sion coefficient, respectively. The subscripts "g" and "d" mean the gas and dust particles, 
respectively. The specific intensity is X v and the mean intensity is J v (= (An)' 1 J X v dQ). 
The first term of the right hand side means the extinction by absorption and scattering, the 
second and third terms are the source term due to the dust thermal continuum emission and 
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the line emission, and fourth term is the source term due to the scattering. We assume that 
the effects of the absorption and the scattering of line emissions by the gas molecules can be 
included into the emission coefficient j Vig and we denote the net emission coefficient as f . 
Using the notation, the radiation transfer equation can be rewritten as 



dx 



= -(a„,d + cr uA )l v + j u>d + j' + a vA J v . (A2) 



Integrating Eq. (A2) over v, we obtain 

— = -(a + a)l + j d +j' g + aJ, (A3) 

where a = J a vA X v dvj j X u du, a = J a vA I v dv / J T u du, and a = J a vA J v dvj j J v dv are 
the frequency- averaged absorption/scattering coefficients. Here, we replace these frequency- 
avaraged absorption/scattering coefficients to the Planck mean values. We think that it is an 
appropriate approximation because the line center frequency is similar to the peak frequency 
of the dust thermal continuum emission. Finally, we can rewrite Eq.(A3) into a simpler form 

^ = -K + a p )(J + <S), (A4) 

where 

S = jj±A±^g. (A5) 

«p + Op 



B. Escape Probability Method 

In the shock-wave heating model, the optical depth effects are important in determining 
the net cooling rate by line emissions of gas molecules. In order to evaluate the net cooling 
rate, we adopt the fitting formula derived by Neufeld & Kaufman (1993), in which they 
numerically solved the equations of statistical equilibrium for the rotational/vibrational level 
populations using an escape probability method. 

The equations of statistical equilibrium is given as 

fMjrh + Cji) ~ fi + (';,). (Bl) 

j 3 

where /j is the fractional population in state i, Ay and CV,- are the spontaneous radiative 
rate and collisional rate for transitions from state i to state j, and (5^ is an angle-averaged 
escape probability and given by 

^ " TT3V < B2 > 
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where the Sobolev optical depth rs is given by 

TS = ^v!/dz\ ifjBji - fiBi ^ (B3) 

where n(M) is the volume density of coolant molecules, B^ and Bij are the absorption and 
stimulated emission coefficients, and dv z /dz is the local velocity gradient. Eq. (B3) is exact 
in the limits of small and large rs, and deviates from the exact expression of Hummer & 
Rybicki (1982) by at most 15%. Solving the equations of statistical equilibrium is made more 
difficult by the fact that the escape probabilities for the emitted line photons are themselves 
a function of the level populations. Neufeld & Kaufman numerically solved those equations 
including a large number of rotational states and obtained the net cooling rates over a wide 
range of physical conditions. 

Neufeld & Kaufman (1993) found a convenient analytic fit to express their numerical 
solutions. The net cooling rate A may be conveniently expressed as a rate coefficient, L, 
defined such that the cooling power per unit volume due to collisions of H2 and species M is 
given by A = Ln(H 2 )n(M), where n(H 2 ) and n(M) are the H 2 and coolant particle densities. 
The cooling rate coefficient L has units of erg cm 3 s _1 . The analytic fit for L is written as 



I = J_ + n ( H2 ) + J_ P( H 2) l a (l_ n V2^0 

L Lq Lute Lq L 711/2 -I V Lute 



(1-=^), (B4) 



where L is a function of temperature alone and Lute, ^1/2, o. are functions of both tem- 
perature and iV(M), and JV(M) is the optical depth parameter. In a case that the emitting 
point is located at the center of a static plane-parallel slab of thickness d, N(M) is expressed 
as 

m - (B5> 

where Av is the velocity dispersion of the coolant molecules. For vibrational cooling, Neufeld 
& Kaufman recommended a two-parameter fit to the net cooling rate which includes only 
the first two terms in Eq. (B4). The values of L , Lute, ^1/2, and a are listed in Tables 2-5 
of Neufeld & Kaufman (1993). 



B.l. Net Cooling Rate 

We consider line photons emitted at x = xa in the post-shock region (see Fig. 15). 
We assume that the photons which achieve at x = x m without being absorbed can remove 
the gas thermal energy to outside. Similarly, it is assumed that the gas molecules in the 
pre-shock region do not absorb the line photons emitted in the post-shock region because 
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of the doppler shift; the relative velocity between the pre- and post-shock gases is larger 
than the velocity dispersion of gas molecules. We define the optical depth parameter of the 
upstream (downstream) side of the emitting point, N up (N down ^ as 

N up (M) = 2 r° mit ^j^-dx, N down (M) = 2 H ^j^-dx. (B6) 

We obtain the rate coefficient for each direction, L up and L down , by substituting N up and 
jydown - m ^ (B4). Finally, we obtain the cooling rate for each direction as 

A u P = L ^ n (H 2 )n(M), A down = L down n(R 2 )n(M). (B7) 

The net cooling rate is given as 

A' = . (B8) 



B.2. Emission Coefficient 

In general, assuming an isotropic emitter, we can write the emission coefficient for the 
radiation transfer equation (Eq. A5) as j g = A/(4n). So, in this study, we obtain the 
emission coefficient as 

(A up /(Att), forx A >x B , 
h \A down /(47r), forx A <x B , 

where A up and A down are the net cooling rates for the direction of upstream and downstream, 
respectively (see Eq. B7). 

However, it would be an underestimation if we give the emission coefficient as A'/ (47r) 
because we assume that only the photons which can escape from the post-shock region 
contribute to the net cooling rate for the gas (see Eq. B8). Let's obtain the mean intensity 
J and net flux T at x = x B in Fig. 15. The net cooling rate A' at x = x A does not include 
the line emissions which are absorbed by gas molecules located from x = to x = x A - On the 
contrary, when we calculate the radiation field at x = x B , we should consider the absorption 
between x = x B and x = x A . Therefore, the mean intensity obtained by our model might 
be weaker than that expected in a real situation. 

How much different is the mean intensity J between our model and the real situation? 
We think that the difference is not significant in the dust temperatures just in front of the 
shock front. Assuming the isotropic radiation field in the pre-shock region and the optically 
thick environment, the dust temperature just in front of the shock front T s f is given by Eq. 
(22), where / indicates the net flux of the radiation field returning from the post- to pre-shock 
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region, in the unit of the gas energy flux flowing into the shock front. For example, when 
all the gas energy flux is converted into the radiation energy at the post-shock region and it 
returns toward the pre-shock region, we obtain / = 1. Therefore, in Fig. 8, a curve labeled 
by "/ = 1" gives the upper limit of T s f under the assumption of the isotropic radiation field. 
Although the dust temperature in a case of the dust model L01 exceeds the upper limit, it 
would be due to an anisotropic radiation field resulting from the optically thin environment. 
In an optically thick environment (r > 1), since the radiation field tends to be isotropic, 
the dust temperature just in front of the shock front would be equal or less than the upper 
limit. The upper limit indicates that when the pre-shock optical depth r pre is less than about 
unity, T sf becomes less than 1273 K and therefore the isotopic fractionation would not take 
place. To summarize, we might underestimate the mean intensity of the radiation field in our 
model, however, our conclusion that the critical optical depth for the isotopic fractionation 
is about 1 — 2 does not change significantly. 

However, there is a possibility that the anisotropy of the pre-shock radiation field affects 
the dust thermal histories in the pre-shock region slightly. Therefore, in order to discuss the 
dust thermal histories with higher accuracy, we need to develop a new radiation transfer 
code for the line emissions. We would like to solve the problem in the future. 

C. Lyman a Cooling 

Regarding the Lya cooling, we neglect the absorption by other hydrogen atoms because 
the optical depth for the Lja photon is small. The absorption coefficient for the Lya emission 
is given as a u = ^-n\B\i${y)i where h is the Planck constant, n\ is the number density of 
hydrogen atom in ground state (n = 1), B i2 is the Einstein B-constant, and <f>{v) is the 
line profile function (Rybicki & Lightman 1979). Einstein B-constant is given from the 
Einstein A-constant by Einstein relation. The line profile function has the maximum value 
of ~ 1/ Au at the line center, where Au is the Doppler broadening in frequency. In the shock 
condition nllv55, n x ~ 10 9 cm -3 and Au ~ 10~ A u behind the shock front, so we obtain 
a v ~ 1CT 5 cm -1 . It indicates that Lya photons are not absorbed until they achieve at about 
10 5 cm behind the shock front. However, our calculation results show that the gas cools until 
it reaches that point for the shock condition nllv55. Moreover, the gas decelerated as it 
cools in the post-shock region. The relative velocity between the shocked gas immediately 
behind the shock front and the cooled gas far from the shock front is about 10 6 cms" 1 and 
it is comparable to the velocity dispersion of the gas molecules. Therefore, the large relative 
velocity does not allow the cooled gas to absorb Lya photons emitted immediately behind 
the shock front. 
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The cooling rate due to Lya emission is given by Spitzer (1978) as 

A Lya = 7.3 x 10" 19 n e "- H exp(-118,400/T)ergcm- 3 s" 1 , (CI) 

where n e and nn are number densities of electron and hydrogen atom in unit cm -3 , respec- 
tively, and T is the gas temperature in unit K. 

D. Constancy of Radiation Flux in Pre-shock Region 

D.l. Reason for Constancy 

The radiation flux JF changing with the distance from the shock front indicates that 
there are some energy sources (for dJ-'/dx > 0) or sinks (for dT jdx < 0) of the radiation. 
We can consider following two possibilities as the radiation sources/sinks; (1) the medium 
(gas and/or dust particles) continues to emit or absorb the radiation energy and (2) the 
medium absorbs the radiation energy and then emits it at other places. 

Let us consider the first possibility. In the post-shock region, it can occur. Behind the 
shock front, there are the hot gas and dust particles which are heated by a shock wave. 
They emit the radiation, then cool and flow downstream, however, the new medium which 
is heated at the shock front continues to be supplied. Therefore, there are eternal radiation 
sources behind the shock front (for the steady state). On the contrary, in the pre-shock 
region, the gas and the dust particles cannot be the radiation energy sources. In this region, 
the main energy source is the radiation coming from the post-shock region. The gas and 
the dust particles in the pre-shock region just emit the radiation energy which is equal to 
the radiation energy that they absorbed. It indicates that they can neither produce nor 
consume the radiation energy inside themselves. Therefore, in the pre-shock region, the first 
possibility will vanish. 

The second possibility will also vanish if we consider the non-relativistic shock velocity. 
The radiation energy coming from the post-shock region diffuses into the pre-shock region 
at the effective speed of ~ c/3r prc , where c is the light velocity in the vacuum and r pre is 
the optical depth of the pre-shock region (e.g., Mihalas & Mihalas 1999, p. 353). Considering 
Tpre ~ 100, the light diffusion velocity is estimated as c/3r ~ 10 3 kms _1 . Since this is much 
larger than the shock velocity of about 7kms~ 1 adopted by Desch & Connolly (2002) as 
the canonical case, the radiation energy can escape from the pre-shock region so fast that 
the gas and the dust particles cannot transport the radiation energy. Therefore, the second 
possibility will also vanish. 

To summarize, the radiation flux in the pre-shock region should be constant. It indicates 
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that the radiative heating takes place even far from the shock front. 

D.2. Reason for Non-Constancy 

In Fig. 5, the radiative flux T of model L10 does not seem to be constant in the pre- 
shock region. We thought that the reason is the shortage of the iteration to converge the 
mean intensity in the computational domain. In our model, we solve the gas/dust dynamics 
and the radiation transfer by turns until the mean intensity does not change more than 1% 
(see sec. 2.1). However, this condition does not guarantee the complete convergence. 

Consider that the number of iteration which is required to satisfy above converge condi- 
tion, iVite, exceeds 100. It indicates that the convergence of the mean intensity is very slow. 
If the iteration process has been continued more and more, the mean intensity would change 
near 1% per iteration for a while. Repeating the iteration further 100 times, the mean in- 
tensity might change ~ 1% x 100 ~ 100%! Therefore, roughly speaking, the above converge 
condition is not enough for the complete convergence when N- lte > 100. The numerical results 
for models of L03 and L01 would have converged enough because N- lte < 100. However, for 
model of L10, the numerical result would not have converged completely because N- ltc > 100. 

However, we do not think that we need more iteration in order to converge completely. 
In model of L10, the isotopic fractionation takes place because the dust temperature just in 
front of the shock front T s f exceeds 1273 K. If we have carried out the more iteration, T s f 
would increase more and more. This does not change the result that the isotopic fraction- 
ation occurs. Therefore, we think that our converge condition is meaningful to obtain our 
conclusions. 

Regarding the radiative flux in the post-shock region, there is another region. In the 
middle panel of Fig. 4, we found that the gas temperature decreases near the right boundary 
of the computational domain (x ~ x m ). From the bottom panel of Fig. 4, we also found 
that the cooling source is the dust particles, namely, the gas thermal energy is converted into 
the dust thermal energy by the thermal collision between the gas molecules and the dust 
surface. Since the dust particles are radiative equilibrium far behind the shock, the thermal 
energy transfered from the ambient gas radiates as the thermal radiation. Since the thermal 
emission works as the radiation source, the radiative flux is not constant at this region. 
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Table 1: Notation of the shock condition. The pre-shock gas number density is n and the 
shock velocity is v s . 



notation 


n 


cm 3 




v s [km s 1 




nl4v08 


10 14 


8 


nl4vl0 


10 14 


10 


nl4vl3 


10 14 


13 


nl3vl6 


10 13 


16 


nl3v20 


10 13 


20 


nl3v25 


10 13 


25 


nl2v30 


10 12 


30 


nl2v35 


10 12 


35 


nl2v40 


10 12 


40 


nllv50 


10 11 


50 


nllv55 


10 11 


55 


nllv60 


10 11 


60 



Table 2: Dust models we adopted in this study are listed. The charactor "P" in the dust 
model means power-law size distribution and "L" means lognormal one. The number "10" , 
"03", and "01" stand for the dust-to-gas mass ratio Cd = 0.10, 0.03, and 0.01, respectively. 
C d . 



notation 


size dist. 


c d 


P10 


power-law 


0.10 


P03 


power-law 


0.03 


P01 


power-law 


0.01 


L10 


lognormal 


0.10 


L03 


lognormal 


0.03 


L01 


lognormal 


0.01 
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Table 3: Calculation results of dust thermal histories of a dust particle of the initial radius 
ao = 100 jum are summarized; the optical depth of the pre-shock region r pre , the dust tem- 
perature just in front of the shock front T s f, the peak dust temperature T pea k, the heating 
rate i?heat, the cooling rate R coo i, and the final dust radius ag n . / indicates a ratio of the net 
flux of the radiation passing though the shock front to the gas energy flux lowing into the 
shock front. a(b) = a x 10 fe . 



shock 


dust 


Tpre 


T sf [K] 


Tp Ca k [K] 


i? heat [K/hr] 


i? cool [K/hr] 




/ 


nl4v08 


P10 


15.7 


1469 


1576 


157 


21.0 


89 


0.459 




P03 


4.58 


1241 


1429 


— 


— 


100 


0.526 




P01 


1.44 


1059 


1329 


— 


— 


100 


0.615 




L10 


2.45 


1153 


1379 


— 


— 


100 


0.616 




L03 


0.73 


1023 


1312 


— 


— 


100 


0.697 




L01 


0.24 


972 


1289 


— 


— 


100 


0.711 


nl4vl0 


P10 


14.1 


1590 





— 


— 





0.395 




P03 


4.82 


1426 


1719 


223 


16.6 


83 


0.547 




P01 


1.56 


1244 


1644 


4.78(6) 


7.23(4) 


99 


0.664 




L10 


2.45 


1369 


1709 


326 


1.77(4) 


95 


0.679 




L03 


0.73 


1192 


1636 


4.72(6) 


9.62(4) 


100 


0.730 




L01 


0.24 


1116 


1610 


4.24(6) 


1.21(5) 


100 


0.735 


nl4vl3 


P10 


10.6 


1709 














0.203 




P03 


4.64 


1551 














0.334 




P01 


1.62 


1389 


2079 


329 


7.37 


61 


0.444 




L10 


2.44 


1607 











0.593 




L03 


0.73 


1419 


2101 


196 


6.72 


18 


0.718 




L01 


0.24 


1305 


2071 


2.07(3) 


4.53(4) 


84 


0.709 


nl3vl6 


P10 


1.44 


1024 


1526 


1.98(6) 


— 


100 


0.761 




P03 


0.42 


915 


1500 


1.50(6) 


— 


100 


0.842 




P01 


0.14 


877 


1492 


1.34(6) 


— 


100 


0.865 




L10 


0.24 


945 


1508 


1.67(6) 




100 


0.898 




L03 


0.073 


929 


1504 


1.59(6) 




100 


0.887 




L01 


0.024 


941 


1507 


1.65(6) 




100 


0.846 


nl3v20 


P10 


1.56 


1219 


1790 


8.46(6) 


2.28(4) 


98 


0.769 




P03 


0.45 


1039 


1754 


7.59(6) 


4.07(4) 


100 


0.789 




P01 


0.14 


956 


1741 


7.28(6) 


4.71(4) 


100 


0.788 




L10 


0.24 


1015 


1753 


7.56(6) 


4.37(4) 


100 


0.848 




L03 


0.073 


928 


1738 


7.19(6) 


4.93(4) 


100 


0.753 




L01 


0.024 


833 


1725 


6.88(6) 


5.40(4) 


100 


0.487 


nl3v25 


P10 


1.67 


1439 


1951 


414 


5.50 


55 


0.680 




P03 


0.48 


1272 


1891 


1.01(7) 


7.56(3) 


93 


0.758 




P01 


0.16 


1184 


1866 


9.32(6) 


1.51(4) 


99 


0.735 




L10 


0.24 


1270 


1892 


1.01(7) 


7.96(3) 


94 


0.885 




L03 


0.073 


1189 


1869 


9.39(6) 


1.44(4) 


99 


0.768 




L01 


0.024 


1028 


1835 


8.29(6) 


2.40(4) 


100 


0.249 



-41 - 



Table 3: Continue. 



shock 


dust 


Tpre 


T si [K] 


Tpcak [K] 


i? heat [K/hr] 


R coo i [K/hr] 


afi n \um\ 

nil |_r~ i 


/ 


nl2v30 

11 l^j V (JU 


P10 


0.13 


880 


1361 






100 


0.919 




P03 

-1 v 7 •_) 


040 


871 

Ul 1 


1 358 






1 00 


906 




P01 


0.013 


856 


1355 






100 


0.725 




L10 


0.025 


844 


1355 






100 


0.604 




L03 


0.0074 


808 


1347 






100 


0.271 




L01 


0025 


816 


1348 






100 


0.239 


n 1 2v35 

11 1 — V *_) ' J 


pi n 

1 ±\J 


1 4 

KJ. 11 


1 054 


1 658 


4 1 3(61 


1 25(4") 


1 00 


927 




P03 

-i y >* ) 


042 


1 050 


1 656 


4 1 1 (fil 

1.111 U J 


1 38(4") 

1 .OOI 1 / 


1 00 

1UU 


926 




P01 


0.014 


1049 


1656 


4.11(6) 


1.43(4) 


100 


0.827 




L10 


0.024 


1056 


1660 


4.15(6) 


1.32(4) 


100 


0.875 




L03 


0.0073 


1015 


1650 


3.99(6) 


1.41(4) 


100 


0.391 




L01 


0.0024 


1022 


1651 


4.02(6) 


1.40(4) 


100 


0.329 


nl2v40 


P10 


0.15 


1219 


1958 


1.17(7) 


1.64(4) 


100 


0.934 




P03 


0.045 


1217 


1957 


1.17(7) 


1.57(4) 


100 


0.938 




P01 


0.015 


1220 


1958 


1.17(7) 


1.55(4) 


100 


0.869 




L10 


0.024 


1228 


1960 


1.17(7) 


1.35(4) 


100 


0.924 




L03 


0.0073 


1207 


1955 


1.17(7) 


1.41(4) 


100 


0.575 




L01 


0.0024 


1203 


1954 


1.17(7) 


1.42(4) 


100 


0.396 


n 1 1 v50 

11 1 1 V iJU 


pi n 


01 3 


808 


1443 






1 00 


940 




P03 


0039 


813 


1443 






100 


0.808 




P01 


0.0013 


803 


1441 






100 


0.469 




L10 


0025 


793 


1441 






100 


0.428 




L03 


00074 


805 


1442 






100 


0.406 




LOl 


00025 


819 


1445 






100 


0.402 


n 1 1 v55 

11 IIVlIU 


pi n 

-1 ±\J 


01 3 


879 


1622 


3 45(61 


4 81 (3") 

I.Oll u / 


1 00 


950 




P03 


0.0040 


887 


1622 


3.47(6) 


6.43(3) 


100 


0.857 




P01 


0.0013 


878 


1620 


3.45(6) 


6.84(3) 


100 


0.506 




L10 


0.0025 


870 


1620 


3.43(6) 


6.38(3) 


100 


0.473 




L03 


0.00074 


879 


1621 


3.47(6) 


6.63(3) 


100 


0.426 




LOl 


0.00025 


893 


1623 


3.45(6) 


6.60(3) 


100 


0.420 


nllv60 


P10 


0.013 


948 


1805 


7.47(6) 


6.67(3) 


100 


0.956 




P03 


0.0041 


958 


1805 


7.50(6) 


7.46(3) 


100 


0.874 




P01 


0.0014 


949 


1803 


7.48(6) 


7.74(3) 


100 


0.531 




L10 


0.0024 


965 


1806 


7.52(6) 


7.09(3) 


100 


0.744 




L03 


0.00073 


951 


1804 


7.50(6) 


7.31(3) 


100 


0.442 




LOl 


0.00024 


966 


1806 


7.52(6) 


7.26(3) 


100 


0.434 
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Table 4: Same as Table 3 except for the initial dust radius a = 251 /im. 



shock 


dust 




T, f [Kl 


pcaK L 


Rho*t fK/hrl 


Rrnn] \K/hl] 
^cooi L / J 


Ofin \uiAl\ 


f 


ni4vuo 


Pin 
r 1U 


1 K 7 
10. / 


1 /l^O 

140» 


lOoU 


i c;7 
10 / 


91 Q 
Zl.O 


9/1 1 
Z41 


n a kg 

u.4oy 




rUo 


4.00 


19/11 
1Z41 


1 /I Qfi 
14o0 






9^1 
Z01 


n 

U.0ZO 




rui 


1 AA 


i n^s 

lUOo 


1 ^38 






9^1 
ZOl 


n fii ^ 

U.010 




Tin 


9 A ^ 
Z.40 


1 1 ^9 
IIOZ 


1 ^87 
loo 1 






9^1 
ZOl 


n ri r 

U.010 




t 

JjUO 


n 7^ 

U. ( O 


1 n99 

1UZZ 


1 ^91 
lOZl 






9^1 
ZOl 


n RQ7 
u.oy i 




t ni 


n 9/i 

U.Z4 


0.71 

y / 1 


1 9Q8 

izyo 






9^1 
ZOl 


n 71 1 

U. ( 11 


n!4vlU 


"Din 

r 1U 


1/11 

14.1 


1 £nn 

159U 


1 C01 

18zl 


z /y 


9.5 ( 


i 1 £ 
140 


U.o9o 




rUo 


4.oZ 


1 /I 9£ 
14Z0 


I 798 

I I Zo 


99Q 
ZZo 


1 £J 7 
10. ( 


9Q/1 
Zo4 


n £/i7 
U.04 1 




P01 


1.56 


1243 


1655 


2 03(6) 


2.88(4) 


250 


0.664 




L10 


2.45 


1368 


1720 


326 


8.49(3) 


247 


0.679 




L03 


0.73 


1191 


1647 


2.01(6) 


3.98(4) 


251 


0.730 




L01 


0.24 


1115 


1622 


1.82(6) 


5.03(4) 


251 


0.735 


nl4vl3 


P10 


10.6 


1709 








o 


0.203 




P0.3 


4 64 


1 551 


2122 


249 


7 49 


70 


D 334 




P01 


1.62 


1389 


2086 


330 


7.30 


211 


0.444 




L10 


2.44 


1606 











0.593 




L03 


0.73 


1418 


2107 


196 


6.74 


167 


0.718 




L01 


0.24 


1304 


2077 


3.00(3) 


2.12(4) 


234 


0.709 


niovio 


pi n 

r 1U 


1 /i/i 

1.44 


i n9/i 

1UZ4 


iooy 


q 9nf 
y .zu^o ) 




9M 
ZOl 


n 7fii 

U. 1 01 




pnq 
rUo 


n /I9 

U.4Z 


qi c; 
y 10 


1^19 
101Z 


7 "iCif ^\ 
1 .OU^O ) 




9M 
ZOl 


n 8/19 

U.o4Z 




pm 

rui 


n 1 A 

U. 14 


877 


i ^n/L 

10U4 


O. i U^OJ 




9^1 
ZOl 


n 8fi^ 
u.ooo 




Tin 


n 9/i 

U.Z4 


Q/l/1 
y44 


1 ^9n 

10ZU 


7 Q9f 

/ .yz^o j 




9^1 
ZOl 


n 80s 
u.oyo 




t n^ 

JjUO 


n n7^ 

U.U 1 O 


098 

yzo 


i ^i fi 

1010 


7 f!1 

1 .Ol^OJ 




9^1 
ZOl 


n 887 
U.OO 1 




t m 


n n9/L 

U.UZ4 


y4u 


1^18 
101O 






9^1 
ZOl 


U.O40 


nl ^rOn 

niovzu 


pi n 

r 1U 


1 ^fi 
1.00 


1 91 Q 

iziy 


1 70S 

i i yo 


O.44^0 J 


8 1 8/^ 
o. lo^o ) 


9/1 Q 

Z4y 


n 7fiQ 
u. i oy 




rUo 


n a ^ 

U.40 


i n^o 
iuoy 


I 7f!1 

I I 01 


o.uy 1^0 ) 


1.0o^4 J 


9^1 
ZOl 


n 78Q 

u. i oy 




pm 

rui 


n 1 A 
u. 14 


you 


1 7/L8 


z.yu^u j 


1 8/1/ 


9^1 
ZOl 


n 788 
U. ( oo 




L10 


0.24 


1015 


1760 


3.09(6) 


1.74(4) 


251 


0.848 




L03 


0.073 


928 


1744 


2.93(6) 


1.95(4) 


251 


0.753 




L01 


0.024 


833 


1731 


2.79(6) 


2.13(4) 


251 


0.487 


nl3v25 


P10 


1.67 


1439 


1981 


414 


5.61 


206 


0.680 




P03 


0.48 


1271 


1923 


4.05(6) 


3.26(3) 


244 


0.758 




P01 


0.16 


1183 


1899 


3.77(6) 


6.27(3) 


250 


0.735 




L10 


0.24 


1269 


1924 


4.07(6) 


3.76(3) 


245 


0.885 




L03 


0.073 


1187 


1902 


3.80(6) 


6.23(3) 


250 


0.768 




L01 


0.024 


1024 


1870 


3.42(6) 


1.01(4) 


251 


0.249 
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Table 4: Continue. 



shock 


dust 


Tpre 


T si [K] 


Tpcak [K] 


i? heat [K/hr] 


R coo i [K/hr] 


afi n \um\ 

nil L( 1 


/ 


nl2v30 

11 l^j V (JU 


P10 


0.13 


871 


1377 






251 


0.919 




P03 

-1 v 7 •_) 


D D4D 


862 


1 372 






251 


906 




P01 


0.013 


848 


1369 






251 


0.725 




L10 


0.025 


834 


1370 






251 


0.604 




L03 


0.0074 


793 


1362 






251 


0.271 




L01 


0025 


801 


1363 






251 


0.239 


n 1 2v35 

11 1 — V *_) ' J 


pi n 


1 4 

KJ. 11 


1 031 

i y >* ) i_ 


1 670 

1U 1 u 


1 83(6"l 


6 82(3") 


251 


927 




P03 

-i y >* ) 


042 


1 027 

1UZ ( 


1 667 

1 Oil 1 


1 82(6"l 


5 26(3") 


251 


926 




P01 


0.014 


1026 


1667 


1.82(6) 


5.67(3) 


251 


0.827 




L10 


0.024 


1032 


1673 


1.84(6) 


5.32(3) 


251 


0.875 




L03 


0.0073 


987 


1663 


1.77(6) 


5.75(3) 


251 


0.391 




L01 


0.0024 


994 


1665 


1.78(6) 


5.71(3) 


251 


0.329 


nl2v40 


P10 


0.15 


1186 


1969 


5.07(6) 


1.22(4) 


251 


0.934 




P03 


0.045 


1183 


1968 


5.07(6) 


6.65(3) 


251 


0.938 




P01 


0.015 


1187 


1968 


5.07(6) 


6.42(3) 


251 


0.869 




L10 


0.024 


1195 


1973 


5.08(6) 


5.91(3) 


251 


0.924 




L03 


0.0073 


1171 


1968 


5.06(6) 


5.94(3) 


251 


0.575 




L01 


0.0024 


1167 


1967 
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Table 5: Same as Table 3 except for the initial dust radius a = 398 /xm. 
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Table 5: Continue. 
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Fig. 1. — A schematic view of our calculation region. Gas and precursor dust particles go 
together at the same speed from pre-shock boundary (x = —x m ). They pass though the shock 
front (x — 0), interact each other in the post-shock region, and finally reach the post-shock 
boundary (x = x m ) and run away to the outer region in the nebula. The radiation emitted 
by the gas and dust particles in the post-shock region is absorbed by the dust particles in 
the pre-shock region. The dust particles in the pre-shock region also absorb the radiation 
emitted by other dust particles in the pre-shock region. In this paper, we call this effect the 
blanket effect. 
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10 100 
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Fig. 2. — Shock conditions for which we perform simulations are displayed as a function 
of the shock velocity v s and the pre-shock gas number density n®. Gray region indicates 
the shock condition which is sufficient to heat precursor dust particles up to 1473 — 2000 K 
(INSN). 
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Diameter [Jim] 

Fig. 3. — Initial dust size distribution model we adopted in this study is plotted in case 
of lognormal one. The vertical axis is normalized cumulative number of chondrules whose 
radii are smaller than that of horizontal axis. The thick solid curve is the lognormal size 
distribution that we adopt in this study. Measured chondrule size distributions are also 
displayed for comparison; each symbol means results of LL3 chondrites by Nelson & Rubin 
(2002) (LL3), L3 chondrites by Rubin & Keil (1984) (L3), EH3 chondrites by Rubin & 
Grossman (1987) (EH3), and C03 chondrites by Rubin (1989) (C03), respectively. 
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Fig. 4. — Abundance of hydrogen molecule (j/h 2 )> g as temperature (T g ), and heating/cooling 
rates for gas in the post-shock region are plotted as a function of distance from the shock 
front. The shock condition is nl4vl0 and the dust model is L01. The unit of the heat- 
ing/cooling rates is ergcm~ 3 s _1 . 



-50- 




Fig. 5. — Radiation fields in the chondrule forming-region for various dust models are plotted 
as a function of distance from the shock front. Shock condition is nl4vl0. Top panel is the 
mean intensity in unit of 10 7 erg cm -2 ster -1 s _1 and bottom panel is the radiation flux in unit 
of 10 7 ergcm~ 2 s _1 . The negative flux means the net radiation flux flows toward upstream. 
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Fig. 6. — Dust thermal histories in the pre-shock region are plotted as a function of distance 
from shock front. Shock condition is nl4vl0, and dust models are L10, L03, and L01. A 
horizontal dotted line indicates 1273 K, above which the isotopic fractionation takes place 
until dust temperature exceeds 1473 K. 
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Fig. 7. — Results of dust thermal histories in the pre-shock region are displayed as a function 
of the optical depth of the pre-shock region: (a) dust temperatures immediately in front 
of the shock front, which are elevated by heating due to the ambient radiation field, and 
(b) the heating rate of chondrules in a temperature range of 1273 — 1473 K, which is the 
region surrounded by two dashed lines in the top panel. Each symbol means different shock 
condition; nl4vl0 (circle), nl3v20 (square), nl2v35 (triangle), and nllv55 (inverse triangle), 
respectively. The isotopic fractionation should occur if the heating rate is lower than about 
10 4 K/hr, which is indicated by a dashed line in the bottom panel. 
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Fig. 8. — Comparison of the pre-shock dust temperature estimated by using the radiative 
diffusion approximation (solid curve) with numerically calculated results (filled symbols) for 
various shock conditions: nl4vf0 (top left), nl3v20 (top right), nl2v35 (bottom left), and 
nllv55 (bottom right), respectively. Triangles indicate the results for the log-normal dust 
size distribution and circles are for power-law one. The parameter / is the fraction of the 
net flux of the radiation returning from the shock front to the gas energy flux flowing into 
the shock front. 
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Fig. 9. — Upper limit of the optical depth of the pre-shock region above which the isotopic 
fractionation should occur is plotted as a function of /, which is the fraction of the radiation 
flux returning from the shock front to the gas energy flux entering into the shock front. 
Labels indicate the shock condition of nl4vl0, n!3v20, n!2v35, and nllv55, respectively. 
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Fig. 10. — The dust thermal histories in the pre-shock region for different spatial scale: 
x m = 10 6 km and Cd = 0.01 (top), and x m = 10 5 km and Cd = 0.10 (bottom). The shock 
condition is nl4vl0 and the dust size distribution is power-law. 
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Fig. 11. — Heating rates and cooling rates of chondrules for a = 100 /xm. Each symbol 
means different gas number density of the pre-shock region; nl4 (circle), nl3 (square), nl2 
(triangle), and nil (inverse triangle), respectively. 
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Fig. 12. — Same as Fig. 11 except for a = 251 /im. 
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Fig. 13. — Same as Fig. 11 except for qq = 398 /im. 



-58- 



(a) 




vol. % 






dense 


100 


80 


porous 





20 


radius 


1.000 


1.216 


cross 
section 


1.000 


1.480 



porosity = 0% 



porosity = 80% 



50 
50 

1.442 
2.080 





100 
1.710 
2.924 



Fig. 14. — Schematic pictures of precursor dust particles in various structures (fraction of 
the dense/porous region). The total masses of those particles are assumed to be the same. 
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Fig. 15. — Emitting point xa and a point at which we calculate the radiation field xb- 



